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Separation  Criteria  for  Three  Dimensional  Steady  Viscous  Flow  /I 


and  Flow  Behavior  Near  Separation  Line 
Zhang  Hanxin 

(China  Aerodynamic  Research  and  Development  Center) 

Abstract 

The  separation  criteria  for  a  three-dimensional  steady 

viscous  flow  are  given  in  this  paper.  The  flow  behavior  near  the 

separation  line  was  investigated.  It  is  pointed  out  that  for  any 

flow  described  by  NS  equations,  the  separation  line  is  the 

"convergent  asymptotic  line"  of  the  wall  limiting  streamline. 

Moreover, it,  itself,  is  also  the  limiting  streamline.  The  beginning 

* 

and  end  of  the  separation  line  were  also  studied  in  this  work. 

It  is  pointed  out  that  either  an  open  or  closed  type  separation 
lineal  extends  towards  infinity  from  the  beginning  and  ends  at  a 
focal  point  or  a  node.  It  cannot  end  at  a  saddle  point.  ' 
Corresponding  conclusions  were  also  obtained  for  the  attachment 
line.  Finally,  a  method  to  determine  the  position  of  the  line  of 
separation  is  discussed. 


I.  Introduction 


It  is  of  important  significance  to  study  the  separation  of 
flow  of  a  viscous  gas  for  the  design  of  aircrafts,  guided 
missiles  and  hypersonic  flying  vehicles.  One  of  the  important 
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issues  is  to  determine  the  beginning  of  separation  and  the 
position  of  the  boundary  of  the  separation  region.  This  is  the 
so-called  separation  criteria  and  separation  line  characteristics, 
with  regard  to  this  problem,  there  are  many  theories, 
experimental  studies [  3]  an(j  comprehensive  reviews 5] # 

Nevertheless,  just  as  Brown  and  Stewartson and  Williams^ 
stress  ,  an  agreement  has  not  yet  been  reached  in  this  area. 
Eichelbrenner  and  Oudartt®!,  Maskell^l,  and  Wang  Guozhangl3! 
believed  that  the  separation  line  is  an  envelope  of  the  wall 
limiting  streamline.  Legendre [71 ,  Lighthill [2] ,  Hunt^8!,  and 
Tobak  and  Peake^-ll]  believed,  however,  that  the  separation  line 
is  a  "convergent  asymptotic  line"  of  the  neighboring  wall 
limiting  streamline.  Moreover, it,  itself,  is  also  a  limiting 
streamline. 

•# 

In  this  paper,  an  attempt  is  made  to  clarify  this  problem 
based  on  theoretical  considerations.  Through  studying  the 
separation  flow,  the  criteria  for  the  separation  of  three- 
dimensional  steady  viscous  flow  are  established!  On  this  basis, 
the  flow  behavior  near  the  separation  line  is  analyzed.  In 
addition,  we  will  prove  that  the  separation  line  is  a  convergent 
asymptotic  line  of  its  neighboring  limiting  streamline. 

Moreover, it,  itself,  is  also  a  limiting  streamline.  Similar 
conclusions  are  obtained  for  attachment.  In  the  paper,  the 
beginning  and  development  of  the  separation  line,  as  well  as  the 
method  to  determine  its  position  are  also  discussed. 


II.  Criteria  for  Flow  Separation 


When  studying  the  separation  of  a  three-dimensional  steady 
flow  on  a  fixed  wall.  Figure  1  is  often  used  to  explain  the 
geometric  significance  of  the  separation  line  in  the  literature: 
(1)  the  separation  line  is  the  intersect  of  the  separated  flow 
surface  and  the  object  surface.  (2)  Fluids  on  both 
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sides  flow  over  the  object,  just  as  the  separation  flow  surface.  /2 
The  criteria  for  separation  are  established  in  the  following 
based  on  this  sense. 


; ,  x 
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/o  7 


Figure  1.  Separation  Picture  and  Coordinate  System 

1.  separation  flow  surface 

2.  streamline 

3.  limiting  streamline 

4.  separation  line 

5.  limiting  streamline 

6.  streamline 

7.  separation  flow  surface 

8.  separation  line 

9.  limiting  streamline 

10.  streamline 


1.  Conditions  for  Separation  Flow  to  Depart  From  Surface 
Let  us  assume  that  x,  y,  z  is  an  orthogonal  coordinate  where 


x  and  y  axes  are  on  the  surface  of  the  object  and  z-axis  is  in 
the  normal  direction.  The  corresonding  scale  indices  are 
hl=hi (x,y, z) ,  h2=h2(x,y,z)  and  h3=l.  It  is  also  assumed  that  the 
equation  of  the  separation  flow  surface  is:  z=f(x,y),  or 

F(x,y, z) =z-f ( x,y) *0  (2.1) 

Then,  the  unit  normal  vector  of  the  separation  flow  is 


n 


i  a/  .  i  df  ;  _ 

h ,  dx  e »~T7  ~djj~e>+e‘ 


df_ 

dy 


(2.2) 


where  ex,  ey,  e2  are  unit  vectors  of. the  x,  y,  z  axes. 

Furthermore,  let  us  assume  the  velocity  of  the  fluid  is: 

V=uex+vey+wez  (2.3) 

Here,  u,  v,  w  represent  the  velocity  components  in  x,  y,  z 

direction.  Due  to  the  fact  that  V  intersects  with  the  flow 

— » 

surface  everywhere,  therefore,  V.n=0.  Based  on  equation  (2.2) 
and  (2.3) ,  we  get 

v  df 

_L  AL.  _  W~  hi  dy  (2.4) 

K  dx - u - 

On  the  other  hand,  the  angle  9  between  the  normal  direction 
of  the  separation  flow  n  and  ex  is: 


We  get 


df  yy_  i  AL 
dy  )  J  A,  dx 


(2.5)* 


♦Here,  we  only  take  the  positive  sign.  The  positive  x  direction 
points  from  the  attached  side  to  the  separated  side. 


By  substituting  equation  (2.4)  into  it,  we  get 


/  3 


V  df 


(2.6) 


This  is  the  equation  determining  the  angle  between  the  flow 
normal  direction  and  the  object  normal  direction.  Now,  let  the 
separation  line  be  the  y  axis  (See  Figure  1) .  It  is  necessary 
for  (n.ey)Q=0.  The  subscript  "0"  represents  the  values  taken 
on  the  separation  line.  By  using  equation  (2.2) ,  we  get 

(2.7) 


=< 

\  h,  dy  It 


This  is  the  geometry  of  the  separation  flow  surface  at  the 
separation  line. 

Equation  (2.6)  is  valid  at  any  point  on  the  separation  flow 
surface  outside  the  separation  line.  Let  o  be  any  point  on  the 
separation  line.  As  shown  in  Figure  1,  let  OG  be  the  intersect 
of  the  coordinate  surface  xoz  and  the  separation  flow  surface  and 


Q  be  a  point  on  it.  Let  us  assume  that  the  arc  length  of  OQ  is 
As.  Let  us  study  the  limiting  form  of  equation  (2.6)  when  as-*0. 
Under  separation  conditions,  the  flow  surface  suddenly  lifts  up 
from  the  object  at  the  separation  line.  Hence,  (tg0)Q^O.  (tg0)o 
is  obtained  by  using  equation  (2.7)  and  the  left  side  of  equation 
(2.6).  However,  due  to  the  fact  that  there  is  no  slipping  on  the 
object  surface,  the  right  side  of  equation  (2.6)  belongs  to  the 
"0/0"  type.  Thus,  we  get  the  following  by  L'Hospital  method: 

■  (t). 


where 


/_d\  =/_d*_\  jL\ 

V  ds  )  o  \  ds  A  \  d x  dx  dz  )a 


represents  the  derivative  in  OG  direction.  Because  Q  is  on  OG, 
(dz/dx) Q=(h^tg0) q.  Thus,  the  above  equation  becomes 

y  ,  I  /  dw  \  .  .  /  dw  \  f  /  V  /  A*.  \  i  /  ■*  , 


<t8  +  «*.  „  ($•). -  [(£).  +  «.  « *,  (-&).](*  #)' 

"“'tl7  ("*T  A)*1'1  ,g  1  nix  "sf)].} 

Using  the  surface  condition  Uq=v0=Wo=0,  we  get 

\  dx  ),  \  dx  /,  V  dx 

Furthermore,  based  on  the  continuity  equation  we  know  that 


(2.8) 


1 

'dip  h,  u)  ,  dip  ht  -j)  ,  dip  k,  ht)  ]  ( 

l  dz  ) 

0 

ip  h ,  h, 

dx  dy  dz  J 1 

Here,  p  is  the  density  of  the  fluid.  It  is  easy  to  find  out  from 
equation  (2.7)  that  the  numerator  of  the  right  hand  side  of 
equation  (2.8)  is  zero.  Hence,  it  is  required  that  (tg0)Q  should 
not  be  zero.  Hence 


(2.9) 


The  expression  for  (tg0)Q  is  further  studied  in  the 
following.  The  origin  of  the  coordinate  can  be  chosen  at  o 
without  losing  generality  because  the  actual  flow  described  by  NS 
equations  has  no  Goldstein^^]  singularity  at  the  separation 
line.  Therefore,  u,  v,  w  can  be  expressed  by  using  Taylor  series 
expansion  near  0. 


(2.10) 


Here,  the  surface  condition,  equation  (2.9)  and  (dw/dz)0=0  are 
utilized.  In  addition,  f(x,y)  and  h2  can  also  be  expanded  as: 


(2.11) 


*  +(-ir). y  +(-fr) 


(2.12) 


In  the  derivation  of  equation  (2.11),  equation  (2.7)  was  applied. 
By  substituting  equations  (2.10) -(2.12)  into  equation  (2.6),  we 


get  the  following  after  some  re-arrangement 


(t60)„  = 


(  dhv'1 

U** , 

1,  ~r 

( 

f  d‘u  \  * 

\dxdz  )t  z 

(2.13) 


where  z=f(x,y).  Based  on  equation  (2.11)  and  (2.5),  we  get 


(2.14) 


(2.15) 


(2.16) 


By  substituting  equation  (2.15)  into  (2.14)  and  followed  by 
substituting  the  result  together  with  equation  (2.16)  into 
equation  (2.13)  and  then  letting  x,  y,  z  approach  zero,  we  get 

’(-fr).a). + K-&-).  - (*•  &).  ]<*.  •«*>.■+ 

+(^).  «■«*»>!-•  <2,17> 


Because  the  origin  o  is  arbitrarily  chosen,  this  equation  is  /5 

valid  at  any  point  on  the  separation  line.  Thus,  the  equation  to 
determine  tg0g  is  obtained.  For  a  Newtonian  fluid,  based  on  NS 
equations,  Newtonian  laws,  surface  conditions  and  equation  (2.9), 
we  have 
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mmmm 


is  the  normal 


Here/  m  is  the  viscosity  index,  p  is  pressure  Tzx 
component  of  the  surface  friction  stress  with  respect  to  the 
separation  line.  By  introducing 


f(y) 


£<v> 


/  d*u  \  /  1  1  dp  \ 

V  dz1  Jo  _  \  M  h ,  dx  ) 


i-w 

m 

5r).“ 

•).  2( 
■(*.$: 

JL.\  (  dv  ) 

,  A,  /«V  dz  ) 

).  <-=- 

L  . 

dr..  \ 
dx  j. 

(*■* 

2( 

dv  W  ht  \ 

dz  /•  \  A,  J. 

dp  \ 
dx  J. 


the  solution  of  equation  (2.17)  is 


(2.19) 


‘'«9>-TS57 


f(y)dy  + 


(A,  tg  9)ci  Jj 


r 


(2.20) 


where  (h^tgejg  represents  the  initial  value  of  (hitge)0  when  y-0. 
Because  equation  (2.9)  is  used,  equation  (2.20)  indicates  that 
(tge)g  has  a  non-zero  solution  when  equation  (2.9)  is  satisfied. 
This  means  that  the  flow  surface  is  lifted  off  the  surface  of  the 
object. 

2.  Recirculation  Conditions 

In  order  to  express  the  second  characteristic  of  the 
separation  flow  shown  in  Figure  1,  we  investigated  the  flow 
behavior  near  the  separation  line.  Its  streamline  equations  are 


(2.21) 


1  dz 

X  XT' 


w 

II 


ft,  dy  v 
h,  dx  "V 


The  first  equation  represents  the  flow  characteristics  in  a 
cross-section  perpendicular  to  the  separation  line.  The  second 
formula  represents  the  sectional  flow  behavior  parallel  to  the 
xoy  plane.  When  z-*0,  it  is  the  wall  limiting  streamline. 

Let  us  substitute  equation  (2.10)  into  the  first  formula  of 
equation  (2.21).  After  neglecting  higher  order  terms,  we  get 


dz 

dx 


(iF)/+2("^r).x 


(2.22) 


Based  on  the  theory  of  singular  point  for  regular  differential 
equations,  when 


/  6 


f  ) 

\  ds1  A1 

{  dxds  ) 

the  point  o  is  a  saddle  point  in  a  cross-section  perpendicular  to 
the  separation  line.  When  q>0,  the  point  o  is  a  nodal  or  focal 
singular  point.  Furthermore,  due  to  the  fact  that  the  flow  is 
away  from  the  surface  of  the  object  separation,  i.e.,  w>0,  it  is 
easy  to  find  out  from  the  third  formula  of  equation  (2.10)  that  (d 
^w/dz^)j>Q.  Hence,  when  q>0,  ( d^u/d  xdz)  q>0  .  Thus, 


*-[(*• -EW-rarlH'' 
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The  node  or  focus  o  is  unstable,  which  is  obviously  not  the  case 
shown  in  Figure  1.  Therefore,  the  point  o  can  only  be  a  saddle 


point.  Thus, 


(d'w  \  (  d-u 
\  3c1  /A  d.v3r 


(2.23) 


Because  O2w/az2)#>0,  therefore, 


(2.24) 


This  is  the  criterion  for  the  fluid  coming  from  both  sides  of  the 
separation  flow  surface  to  move  outward  with  the  separation  flow. 
We  call  it  the  recirculation  conditions. 

3.  Condition  for  Flow  Separation 

Based  on  the  above  discussion,  the  condition  for  determining 
flow  separation  on  the  separation  line  is 


(2.25) 


Similarly,  the  condition  for  attachment  of  a  separation  flow  in 
the  attachment  line  is 


(2.26) 
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III.  Flow  Behavior  Near  Separation  Line  and  Attachment  Line 


Now  let  us  study  the  third  formula  in  equation  (2.21) .  When 
z-»0,  it  gives  the  flow  behavior.  Based  on  the  separation  and 
attachment  criteria  given  above,  we  have  the  following 
conclusions: 

1.  If  (dvlds) :#o,  on  the  separation  line,  the  separation 
line  is  a  wall  limiting  streamline. 

In  reality,  let  us  assume  that  the  wall  limiting  streamline 
at  a  point  o  on  the  separation  line  is  x=x(y) .  In  the  following, 
we  will  find  various  derivatives  of  x(y)  with  respect  to  y  at 
point  o.  Based  on  the  definition  of  the  wall  limiting  streamline, 
the  first  order  derivative  is 


h i  dx  __  d >4 1 4* 

hi  '7y”"*$v7ai  (3.1) 


This  equation  is  used  to  find  its  directional  derivative  along  /7 
the  limiting  streamline: 


A.  d'x  ,  [  d  (  hi  \  ,  d  (  kx  \  dx  .  dx 

Hi  dy1  l  dy  \  h,  J  dx  \  h,  )  dy  J  dy 

d*u  d‘u  dx  V  dv  \~‘  l  du  V  3^ 

dydx  dzd x  dy  A  dz  )  \  dz  J\  dydz 


(3.2) 


Based  on  equation  (2.9) ,  on  the  separation  line,  3u/3z=0.  Hence, 
32u/3y3z=0,  3x/3y=0.  From  equation  (3.2),  we  get 


d'x 

dy1 


Similarly,  we  can  prove  that 


(3.3) 


dr*_ 

dy‘ 


(3.4) 
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which  indicates  that  all  derivatives  of  the 


Here,  n=2,3,4... 

limiting  streamline  x*x(y)  are  zero  at  point  o.  Therefore,  this 
limiting  streamline  is  the  y  axis;  i.e.,  the  separation  line  is  a 
limiting  streamline. 

It  is  also  possible  to  prove  that  this  conclusion  is  correct 
from  another  angle.  In  reality,  we  can  perform  Taylor  series 
expansion  of  au/az  and  av/az  at  point  o: 


(3.5) 


After  substituting  it  into  equation  (3.1)  and  neglecting  higher 


order  small  terms,  we  get: 


(3.6) 


where 


A  = 


(JLl\  (  d>uldxd:  \ 

\  A,  /.\  dv/dz  /• 


(3.7) 


Integrating  equation  (3.6)  we  get 

jc  —  ±ceJt  (3.8) 

Here,  c  is  a  constant  of  integration.  It  is  possible  to  see  that 
if  the  integration  curve  is  required  to  pass  through  the  point  o 
(x*y=0) ,  c=0.  This  means  that  the  limiting  streamline  passing 
through  o  is  x=0.  Therefore,  the  separation  line  is  a  limiting 
streamline. 

2.  If  av/az=0,  then  the  limiting  streamlines  in  the 
neighborhood  of  the  separation  line  converge  toward  the 


separation  line.  Moreover,  the  separation  line  is  the  asymptotic 
line. 

In  reality,  based  on  Conclusion  1,  equation  (3.8)  is  valid. 

As  shown  in  Figure  2,  let  the  x  coordinate  for  the  point  A  be  e  . 
Then  the  limiting  streamline  passing  through  A  is  x=eeAY.  Near 
the  surface  of  the  object,  v  >0.  Thus,  the  flow  is  in  the  same 
direction  as  the  positive  y  axis  direction.  Because  v=0  on  the 
surface  of  the  object,  9v/3z  >0.  For  the  separation  line, 
however,  32u/9x9z  <0.  Based  on  equation  (3.7)  we  know  that  A  <0. 
This  indicates  that  the  distance  between  the  limiting  streamline 
and  the  separation  line  decreases  with  increasing  y.  Moreover, 
when  s^m.x-O.  This  means  that  the  limiting  streamline  is 
converging  towards  the  positive  y  direction.  If  v<  0,  then  the  /8 
flow  direction  agrees  with  the  negative  y  direction.  In  this 
case,  3v/3z  <0  on  the  surface  of  the  object.  Therefore,  for  the 
separation  line,  the  distance  between  the  limiting  streamline  and 
the  separation  line  decreases  with  decreasing  y  when  A>  0.  In 
addition,  when  y-  x-0.  The  limiting  streamline  converges  in 
the  negative  y  direction.  The  study  shows  that  in  all  cases  the 
separation  line  is  the  converging  line  of  the  neighboring 
limiting  streamline. 

3.  If  3v/azH0,  in  the  neighborhood  of  the  attachment  line, 
the  wall  limiting  streamline  diverges  using  the  attachment  line 
as  the  asymptotic  line  (see  Figure  3). 


Figure  2.  Behavior  of  Limiting  Streamline  Near  Separation  Line 
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Figure  3.  Behavior  of  Limiting  streamline  Near  Attachment  Line 

1.  limiting  streamline 

2.  attachment  line 

3.  attachment  line 

4.  limiting  streamline 


As  shown  in  Figure  3,  the  equation  for  the  limiting 
streamline  passing  through  A  is  x=eeA¥.  If  v  >0,  then  av/az  >0. 
In  the  event  of  attachment,  a^u/axaz  >0.  Therefore,  when  A>  0, 
the  distance  between  the  limiting  streamline  to  the  separation 
line  increases  with  increasing  y.  When  y-  -®,  x-0.  The  limiting 
streamline  diverges  along  the  positive  y  axis.  If  v  <0,  then 
av/az  <0  and  A  <0.  When  y-«,  x-0.  The  distance  between  the 
limiting  streamline  and  the  separation  line  increases  as  y 
decreases.  This  means  that  the  limiting  streamline  diverges 
along  the  negative  y  axis.  In  summary,  regardless  of  the 
situation,  the  limiting  streamline  near  the  attachment  line 
diverges  outward  using  the  attachment  line  as  the  asymptotic 
line. 

4.  For  a  real  flow  satisfying  the  Newtonian  friction  law, 
the  separation  criterion  is: 


The  attachment  criterion  is: 


(3.10) 


Furthermore,  conclusions  1,  2  and  3  are  also  applicable  to  the 
friction  line. 

This  is  due  to  the  fact  that  based  on  Newtonian  law  that  on 
the  surface  of  the  object 


The  friction  line  is  the  limiting  streamline. 


IV.  Beginning  and  Development  of  the  Separation  Line 

In  general,  the  most  upstream  position  of  the  separation 
line  is  the  beginning  of  the  separation  line.  Without  losing 
generality,  through  a  proper  choice  of  the  coordinate,  there  are 
two  situations  concerning  the  beginning  of  the  separation  line. 

One  is  a  normal  point  where  both  friction  components  are  non¬ 
zero.  The  other  is  a  singular  point  where  both  friction  / 9 

components  are  simultaneously  zero.  Based  on  the  theory  that  a 
unique  solution  to  a  first  order  normal  differential  equation 
always  exists,  the  separation  line  starting  from  a  normal  point 
is  the  segment  of  the  limiting  streamline  which  passes  through 
the  starting  point.  The  condition  represented  by  equation  (2.25) 
is  met  on  it.  Furthermore,  (d‘u/dxdz),< o  distinguishes  this  segment 
from  other  neighboring  limiting  streamlines.  This  is  the  open 
type  of  separation [31  described  by  Wang  Guozhang  (see  Figure  4a). 
The  behavior  of  the  separation  line  beginning  from  a  singular 
point  will  be  analyzed  together  with  other  singular  points  on  the 
separation  line. 
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Figure  4.  Behavior  of  Beginning  of  Separation  Lines 


An  orthogonal  coordinate  x,  y  is  chosen  on  the  surface  of 
the  object.  (Here,  the  y-axis  needs' not  to  be  the  separation 
line.)  The  equation  of  the  limiting  streamline  can  be  written 
as: 


h 

dy  1  dz 


dx  ,  du 

hl~dT 


(4.1) 


Assume  that  s  is  a  singular  point  on  the  separation  liner 

Mr). 


«■  -  [£(*•  -£ )  •  £(*■  -=-)-£(*■  -It)  •  £(*•  -I-)] . 
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(4.2) 
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Based  on  the  theory  of  singularity  in  differential  equations, 
when  qs < 0,  it  is  a  smaller  point.  When  qs  > 0,  the  singular 
point  is  a  node  or  a  focus.  Furthermore,  when  Rs  >  0,  it  is 
stable  and  the  limiting  streamline  points  toward  the  singular 
point.  When  Rs<0,  it  is  unstable  and  the  limiting  streamline 
runs  away  from  the  singular  point.  Based  on  these  results,  the 
following  conclusions  can  be  obtained: 

1.  If  the  start  of  the  separation  line  is  a  singular  point, 
it  must  be  a  saddle  point. 

As  a  matter  of  fact,  based  on  the  second  formula  of  equation 
(4.2)  and  the  continuity  equation  we  get 


R, 


(4.3) 


Because  a2w/az2  >0  during  separation,  therefore,  Rs  >0.  This 
indicates  that  if  the  singularity  is  a  mode  or  focus,  it  will  be 
stable.  This  is  obviously  not  the  beginning  of  the  separation 
line  because  the  direction  of  the  limiting  streamline  at  the 
beginning  of  separation  is  outward.  In  order  to  prove  that  the  /10 
singularity  is  indeed  a  saddle  point,  let  us  choose  the  y-axis  as 
the  separation  line.  Because  a u/9 z=0  along  the  y-axis,  therefore 
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In  addition,  because  it  is  located  on  the  separation  line,  we 
have 


d‘u 

dxdz 


<0 


Furthermore,  it  is  at  the  start  of  the  separation  line 


Thus,  qg <  0  which  corresponds  to  a  saddle  point.  In  this  case, 
the  separation  line  starts  from  the  saddle  point  and  propagates 
in  both  directions.  Moreover,  the  two  branches  intersect 
tangentially  at  the  starting  point,  forming  a  smooth  curve 
extending  .on  both  sides  and  separating  the  incoming  flow  from  the 
separation  region  (see  Figure  4b) .  This  is  the  so-called  closed 
type  of  initial  separation  described  by  Wang  Guozhangt3!. 

2.  If  in  the  developing  process  the  separation  line 
encounters  singular  points,  they  can  only  be  nodal  or  focal 
points. 

As  a  matter  of  fact,  as  the  separation  line  enters  a 
singular  point,  we  have 


However,  based  on  equation  (2.25) 


Hence,  based  on  the  first  formula  in  equation  (4.4),  qs>0. 

Thus,  it  is  not  possible  to  be  a  saddle  point. 

Furthermore,  equation  (4.2)  can  be  used  to  analyze  that  the 


singular  point  may  be  a  nodal  or  focal  point. 

Similarly,  the  attachment  line  can  be  studied.  It  was  found 
that  it  also  might  have  two  patterns.  One  starts  from  a  normal 
point.  The  other  begins  from  a  node  or  a  focus.  It  is  not 
possible  to  start  from  a  saddle  point.  The  end  is  a  singular 
point  which  must  be  a  saddle  point. 

Finally,  it  should  be  pointed  out  that  when  z-0,  at  the 
singular  point  s(x=y=0),  the  second  formula  in  equation  (2.21) 
can  be  expressed  as  the  following  by  using  NS  equations: 

1  /  d'y  \  J  /  d‘v  \ 

(fhdy  \  m  TA  dz‘),2  ~t~"‘  _  V  dz:  ),  _(h_x\(  d_pjdy  \ 

•\hxdx),  l/aiu\  ,,  (  d:u\  Vh'JA  dp/dx  },  (4-5) 

2  \  dz1 ) ■ 2  •  Va?-/. 

It  is  possible  to  see  that  the  pressure  p  at  point  s  is  not  an 
extremum;  i.e.,  (ap/ax)  and  (ap/ay)  are  not  simultaneously  zero. 
Equation  (4.5)  gives  a  specific  direction  for  the  curve. 

Obviously,  it  is  not  the  spiral  singularity  observed 
experimentally.  Therefore,  the  necessary  condition  for  the 
spiral  singularity  is  that  pressure  is  at  its  extreme  at  s. 

V.  Determination  of  Separation  Line  Position  /II 

In  the  literature,  the  position  of  the  separation  line  is 
often  determined  by  the  two  following  methods.  One  is  to  find 
the  converging  asymptotic  line  or  envelope  of  the  wall  limiting 
streamline ^ ^ 1 .  Based  on  the  above  analysis,  only  the  former 
has  some  basis  for  flows  described  by  NS  equations*.  The  other 


^  a 


is  based  on  the  minimum  surface  friction  line  and  assumes  that 


the  separation  line  is  the  converging  asymptotic  line  of  the 
limiting  streamline  which  is  parallel  to  the  minimum  friction 
line.  Then,  through  series  expansion  of  the  friction  force 
equation,  the  position  of  the  separation  line  is  determined  by 
taking  the  lower  order  terms  1^*1.  Based  on  the  separation 
criteria  described  in  this  paper,  the  assumption  in  this  method 
can  be  dropped.  Thus,  the  discussion  in  this  paper  provided  a 
theoretical  basis  for  this  method.  Furthermore,  this  method  is 
also  expanded. 

This  method,  in  principle,  is  applicable  to  the 
determination  of  the  position  of  the  attachment  line. 


VI.  Simple  Conclusions 


In  a  real  three-dimensional  viscous  flow  described  by  NS 
equations,  theGoldstein  singularity  no  longer  exists  in  the 
neighborhood  of  separation  and  attachment  lines.  Based  on  this 
characteristic,  through  the  analysis  above,  we  arrived  at  the 
following  conclusions: 

1.  The  criterion  for  flow  separation  is:  on  the  separation 

line 


The  criterion  for  the  attachment  of  a  separated  flow  is:  on  the 
attachment  line 

Mr).  (&> 
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This  is  an  extension  of  the  two-dimensional  Prandtl  criterion. 

2.  The  separation  line  is  a  limiting  streamline.  It  is  the 
converging  asymptotic  line  for  neighboring  limiting  streamlines. 
It  is  not  the  envelop  of  limiting  streamlines. 

3.  The  attachment  line  is  a  limiting  streamline.  It  is  a 
diverging  asymptotic  line  for  its  neighboring  limiting 
streamlines.  The  attachment  line  is  not  the  envelope  of  limiting 
streamlines. 

4.  A  separation  line  may  have  two  starting  modes.  One 
begins  with  a  normal  point.  In  this  case,  the  separation  line  is 
the  segment  of  limiting  streamline  beyond  the  starting  point. 

This  is  equivalent  to  the  open  separation  described  by  Wang 
Guozhangl3!.  The  other  type  starts  from  a  saddle  point.  The 
separation  line  is  the  boundary  between  the  incoming,  flow  region 
and  the  separation  region.  This  is  equivalent  to  the  closed 
separation  mentioned  by  Wang  Guozhang.  In  addition  to  extending 
to  infinity,  a  separation  line  may  end  at  a  node  or  a  focus.  If 
the  beginning  of  a  separation  line  is  a  singular  point,  it  must 
be  a  saddle  point. 

5.  An  attachment  line  may  also  have  two  modes  -  normal 
point  and  singular  point.  However,  it  may  only  begin  as  a  focus 
or  a  node.  In  adding  to  extending  to  infinity,  it  can  only 
terminate  at  a  saddle  point. 

6.  It  is  possible  to  find  the  converging  wall  limiting 
streamline  in  order  to  determine  the  position  of  the  separation 
line.  It  is  also  possible  to  start  from  the  minimum  friction 
line.  The  distance  between  the  separation  line  and  the  friction 


line  can  be  found  based  on  the  separation  criteria  to  determine 
the  position  of  the  separation  line.  Similar  methods  can  be  used 
to  determine  the  position  of  the  attachment  line. 

*In  another  paper  we  already  proved  the  envelope  criterion 
is  correct  for  a  flow  described  by  boundary  layer  equations. 

Professor  Liu  Moji  of  Beijing  Institute  of  Aeronautics  and 
Astronautics  is  acknowledged  for  his  assistance. 
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Aerodynamic  Calculation  of  Airfoil  with  Separation  / 1 3 

Zhu  Ziqiang,  Chen  Bingyong  and  Zhang  Bingxuan 
(Beijing  Institute  of  Aeronautics  and  Astronautics) 

Abstract 

An  aerodynamic  computation  for  wing  section  with  separation 
is  presented  in  this  paper.  A  vortex  is  used  to  simulate  the 
separation  wake  boundary.  The  strength  of  this  vortex  sheet  is 
rationally  chosen.  An  "equivalent  body"  is  formed  by  this 
separated  wake  vortex  sheet  and  the  attached  flow  region  of  the 
wing  section.  The  iterative  potential  flow  solution,  which 
includes  boundary  layer  effect,  is  used  to  solve  this  equivalent 
body.  The  shape  of  the  separation  wake  boundary  and  the  position 
of  the  separation  point  are  determined  in  the  iteration  process. 

Two  classical  wings  were  calculated.  The  computational  results 
are  in  good  agreement  with  experimental  data. 

I.  Introduction 

In  most  aeronautical  problems  involving  high  Reynolds 
numbers,  when  the  angle  of  attack  is  not  too  large  and  flow 
separation  does  not  occur,  the  aerodynamic  characteristics  of  a 
viscous  flow  around  an  object  (such  as  a  wing)  can  be  calculated 
with  sufficient  accuracy  through  a  potential  flow  calculation 
with  boundary  layer  effect.  Reference  [1]  presented  the  problem 
of  calculating  the  mutual  interaction  between  the  potential  flow 


and  the  boundary  layer.  The  computational  results  obtained  in 
that  paper  were  found  to  be  in  good  agreement  with  the 
experimental  data  at  small  attack  angles.  With  increasing  attack 
angles,  the  difference  gradually  increases.  This  is  due  to  local 
separation  on  the  wing.  As  the  attack  angle  gets  larger,  this 
separation  expands  from  the  rear  towards  the  front.  Hence,  in 
order  to  rationally  estimate  the  aerodynamic  characteristics  at 
large  attack  angles,  it  is  necessary  to  take  the  effect  of 
separation  on  the  wing  into  consideration. 

As  the  range  of  an  attack  angle  used  in  aerospace  expands,  also 
in  response  to  the  need  in  designing  high  lift  equipment,  it  is 
desirable  to  have  a  practical  engineering  calculation  method 
which  can  account  for  the  effect  of  separation.  It  should  be 
capable  of  calculating  the  aerodynamic  characteristics  within  the 
attack  angle  range  up  to  near  the  maximum  lift  coefficient  Clmax, 
including  the  capability  to  provide  a  relatively  accurate 
pressure  distribution. 

Jacobi^]  first  addressed  this  problem  by  empirically  placing 
rational  source  and  sink  distributions  in  the  trailing  edge  of 
the  wing  to  simulate  wake  separation.  Then,  he  employed  the 
potential  flow  theory  to  solve  the  flow  of  this  "equivalent  body" 
with  distributed  sources.  Results  indicate  that  it  is  feasible 
to  calculate  the  effect  of  separation  wake  by  extending  the 
streamline  in  the  rear  of  the  wing  using  source  and  sink. 
Masanori,  Hayashi  and  Eiichi  Endo^I  also  used  a  similar  method 
to  calculate  the  Clmax  characteristics. 

It  is  also  possible  to  use  a  vortex  sheet  extending  from  the 
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wing  as  the  separation  wake  boundary,  instead  of  using  a  source 


and  sink  distribution.  An  "equivalent  body"  is  formed  by  adding 


such  a  vortex  sheet  in  front  of  the  attachment  flow  region  of  the 


wing.  The  boundary  of  the  separation  wake  and  the  position  of 


the  separation  point  are  determined  in  the  iterative  process  in 
finding  the  complete  solution l4-"?]  .  in  addition,  other  authors 


took  a  step  further  to  consider  the  flow  in  the  separation  wake. 


An  extended  jet  mixing  model  was  used  as  the  solution  of  the  flow 


inside  the  wake.  Then,  it  is  used  simultaneously  with  the 


external  equivalent  potential  flow  solution  to  determine  the 


boundary  of  the  separation  wake.  In  this  case,  iterations  are 


still  required.  Consequently,  another  iterative  step  is  added  to 


complicate  the  entire  process.  The  work  load  is  significantly 
increased^6-?] . 


In  this  paper  is  an  iterative  method  for  calculating  the 


aerodynamic  force  by  using  a  vortex  sheet  to  simulate  the 


separated  wake  boundary.  It  is  capable  of  taking  the  effect  of 


separation  and  the  wall  boundary  effect  at  the  front  of  the  wing 


into  account.  The  intensity  of  the  separation  vortex  sheet,  to 


some  extent,  can  be  rationally  selected  based 
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on  the  pressure  distribution  in  the  separation  wake  area  obtained  / 1 4 


experimentally.  Thus,  the  computational  results  would  be  in 


better  agreement  with  the  reality. 


The  author  wishes  to  express  his  gratitude  to  Professor  Lu 


Shijia  for  his  concern  and  support. 


II.  Aerodynamic  Simulation  of  Flow  Around  Wing  with  Separation 

Wake 

Based  on  a  great  deal  of  experimental  work  on  different 
wings  at  various  attack  angles,  in  the  separation  wake  region, 
the  pressure  distribution  in  front  of  the  trailing  edge  of  the 
wing  (such  as  in  area  A  in  Figure  1)  is  essentially  constant.  It 
is  equal  to  the  pressure  at  the  separation  point.  Afterward, 
during  the  merge  of  the  two  flows,  pressure  will  gradually 
increase  to  CpR  at  the  merging  point  R.  CpR  is  slightly  less 
than  the  pressure  at  the  stationary  point.  We  choose  the  vortex 
sheet  as  the  boundary  of  the  separation  wake  region.  In 
addition,  based  on  the  fact  that  this  vortex  sheet  is  a 
streamline,  this  streamline,  together  with  the  wing  surface  in 
front  of  the  separation  point  where  the  flow  is  attached,  is 
considered  as  an  equivalent  body.  From  the  potential  flow 
solution  of  the  equivalent  body  which  includes  boundary  layer 
effect,  it  is  possible  to  obtain  the  pressure  distribution  on  the 
wing.  Due  to  the  fact  that  the  separation  vortex  sheet  on  the 
upper  wing  surface  begins  at  the  separation  point  of  the  boundary 
layer,  the  profile  of  the  enclosed  wake  boundary  and  the  position 
of  the  separation  point  must  be  solved  by  iterations  as  a  part  of 
the  whole  solution. 

When  there  is  a  potential  flow  around  the  equivalent  object, 
there  should  be  two  stationary  points.  We  choose  the  merge  point 
as  the  rear  stationary  point.  This  is  equivalent  to  the 
situation  that  two  vortex  streamlines  from  the  upper  wing  surface 
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and  the  trailing  edge  meet  at  the  merge  point;  where  the  velocity 


should  be  zero  and  the  corresponding  vortex  strength  should  be 


zero.  Of  course/  the  pressure  is  the  stationary  point  pressure. 


Compared  to  the  pressure  at  the  point  in  a  real  flow,  there  is 


some  difference,  yet  the  difference  is  very  small.  The 


difference  will  not  significantly  affect  the  pressure 


distribution  on  the  wing.  In  region  A  of  the  separation  wake,  we 


choose  the  pressure  at  the  upper  wing  separation  point  and  that 


at  the  lower  wing  trailing  edge  to  be  equal  to  the  pressure. 


When  solving  the  potential  flow,  we  distribute  vortices  on  the 


attachment  flow  wing  surface  and  the  separation  boundary. 


Because  under  the  condition  that  the  object  surface  cannot  be 


penetrated,  the  vortex  intensity  distribution  7  on  an  enclosed 


vortex  sheet  without  any  singular  points  is  equal  to  the  velocity 


outside  the  surface.  The  condition  that  the  pressure  at  the 


separation  point  on  the  upper  surface  is  equal  to  that  at  the 


trailing  edge  of  the  lower  surface  is  equivalent  to  the  fact  that 


the  vortex  intensities  at  these  two  places  are  equal. 


Furthermore,  by  referring  to  the  direction  of  flow,  it  is 


possible  to  see  that  the  vortex  strengths  are  equal  and  the 


directions  are  opposite  at  these  two  places;  i.e.,  7. 


separation=-7low,  trailing  edge*  This  may  also  be  called  the 
extended  Kuttar  condition  with  separation. 


Thus,  the  strength  of  the  separation  vortex  sheet  can  be 


chosen  as  follows.  In  region  A,  the  vortex  sheet  strength  is 


equal  to  that  at  the  separation  point.  In  region  B,  it  is 


decreasing  with  the  square  of  distance  downstream  until  7=0  at 


w 


vVvVv'/  •' 


the  point  of  merge  (See  Figure  1) .  At  the  upper  boundary,  7=-y 
l [l-(x^/xw) 2] .  Here,  7^  is  the  vortex  strength  at  the  lower 
trailing  edge.  The  shape  of  the  separation  vortex  sheet, 
however,  should  be  determined  based  on  its  streamline  condition 
in  iterations. 

As  for  the  calculation  of  the  boundary  layer  displacement 
thickness  effect  of  the  attachment  flow,  similar  to  what  was  done 
in  reference  [1],  the  equivalent  normal  penetration  velocity  can 
be  determined  based  on  the  boundary  layer  displacement  thickness 
derived  from  the  potential  flow.  It  is  only  necessary  to  deploy 
a  corresponding  source  distribution  with  strength  m=(d/dx)  (ue<5*) 
to  simulate  the  boundary  layer  effect. 

III.  Calculation  of  Potential  Flow 

In  order  to  facilitate  the  calculation,  the  upper  and  lower 
wing  surfaces  are  divided  into  many  equal  and  unequal  elements. 
The  polygon  formed  by  these  surface  elements  is  used  to 
approximate  the  profile  of  the  original  object.  Considering  that 
the  aerodynamic  characteristics  near  the  trailing  edge  vary  more 
rapidly  at  low  flow  velocities,  the  following  unequal  method  is 
used  to  divide  the  surface  elements: 

*./c=-i~Cl+coS5,]  (1) 

0,=»'jr/ni  i  =  0, 1,2, ,n 

where  the  origin  is  located  at  the  leading  edge,  x-axis  is  along 
the  chord  direction,  z-axis  is  positive  upward,  and  c  is  the 
chord  length. 


On  small  surface  elements,  linear  strength  vortex 
distribution  is  arranged.  Its  strength  is  determined  by  the 


boundary  condition  that  no  penetration  is  allowed  at  the  control 
point  of  each  element. 

The  induced  velocity  produced  by  any  vortex  segment  B  at  any 
point  in  space  can  be  calculated  based  on  Biot-Savart  equation: 


(2) 


The  normal  induced  velocity  at  point  i  due  to  the  vortex 
distribution  of  the  entire  wing  is: 


(3) 


where  Anifj  is  the  aerodynamic  force  effect  coefficient. 

If  we  choose  a  typical  surface  element  j  and  control  point  i 
(such  as  shown  in  Figure  2) ,  the  vortex  strength  on  the  element 
is  a  linear  function,  i.e., 


•  l 


(4) 
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Figure  1 


Figure  2 


After  substituting  it  into  the  Biot-Savart  equation,  we  get 


Zv.iH,,.,? -h>, 7  7  7 >3  (5^ 

where  H-i  ,  ,  Ho.  ,  Ho.  and  .  .  ,  were  defined  in  reference 
1i  •  J  *i’j  Ji’j  *i-j 

[1]  • 


Comparing  equation  (5)  to  equation  (3) ,  we  can  write 


Ai-i-H f  -  H  7  +  i  "^>,.,7 


In  addition,  the  normal  velocity  component  of  the  incoming  flow 
at  point  i  is: 


P'.  =F»[cos  a  i  -“-sin  a  7  3  * 


Tm*1 


In  order  to  satisfy  the  boundary  condition,  we  get  the  following 
matrix  equation: 


[A]{y}  =  - 

2-7 


(8) 


Based  on  this,  it  is  possible  to  obtain  the  unknown  vortex 
strength  distribution: 

{v}=  -2nlA.y'{V,J  (9) 

From  equation  (8)  we  know  that  the  number  of  control  points 
is  one  less  than  that  of  the  unknown  y.  In  order  to  close  the 
equations,  the  Kuttar  condition  can  be  used  as  a  supplemental 
equation.  In  reference  [1],  several  commonly  used  Kuttar 
conditions  were  compared.  In  this  paper,  because  the  effect  of 
separation  wake  is  considered,  we  adopt  the  following: 

Kip,  separation=“^low,  trailing  edge 
IV.  Iterative  Calculation  to  Consider  Separation  Wake  Effect 

Based  on  the  aerodynamic  model  of  separation  wake  discussed 
above,  the  specific  steps  for  computation  are: 

(1)  Make  potential  flow  calculation  for  the  wing  first  to 
obtain  the  potential  flow  pressure  distribution. 

(2)  Calculate  the  boundary  layer  based  on  the  pressure 
distribution  obtained.  In  addition,  the  displacement  thickness 
distribution  and  the  separation  point  position  on  the  wing  are 


also  obtained. 

(3)  Assume  an  initial  profile  for  the  wake  boundary. 

(4)  Adjust  the  wake  boundary  leaving  the  separation  point 
unchanged.  Based  on  the  requirements  that  the  wake  boundary  is  a 
streamline  and  the  solid  wall  attachment  flow  area  is  not 
penetrable  (including  the  displacement  thickness  effect  of  the 
boundary  layer) ,  iterations  are  performed  until  the  wake  boundary 
variation  becomes  very  small  after  two  iterations.  We  call  it 
internal  iteration. 

(5)  The  new  pressure  distribution  thus  obtained  is 
different  from  the  originally  calculated  pressure  distribution  in 
the  boundary  layer.  It  is  necessary  to  re-calculate  the  boundary 
layer  based  on  the  new  pressure  distribution  in  order  to  obtain 
the  new  boundary  layer  displacement  thickness  distribution  and 

* 

the  separation  point  position. 

(6)  Repeat  steps  (3) -(5)  as  the  iterative  calculation.  We 
call  this  external  iteration.  It  is  carried  out  until  the  change 
of  in  two  consecutive  iterations  is  less  than  a  specified 
number.  The  convergent  solution  at  this  attack  angle  is  thus 
obtained. 

In  the  calculation,  because  the  separation  point  position 
varies,  it  is  necessary  to  make  geometric  modification  for  each 
new  separation  point  position  in  order  to  make  the  separation 
point  locate  at  a  corner  of  the  surface  element.  Hence,  we  must 


also  modify  a  portion  of  the  aerodynamic  effect  index  matrix.  In 
addition,  the  changing  vortex  sheet  also  requires  the 
modification  of  the  aerodynamic  index  matrix  accordingly. 

However,  only  individual  rows  or  columns  are  required  to  be 
modified.  It  is  not  necessary  to  repeat  the  calculation  of  the 
entire  aerodynamic  index  matrix  in  order  to  save  computing  time. 
Thus,  in  every  step  of  the  iterative  process,  after  modifying  the 
aerodynamic  effect  index  matrix,  new  vortex  strength  distribution 
and  pressure  distribution  that  satisfy  the  boundary  condition  of 
the  wing  surface  can  be  obtained. 

In  the  iteration,  the  third  step  is  to  provide  an  initial 
profile  for  the  wake  boundary.  In  this  calculation,  it  is  given 
according  to  the  following.  The  wake  boundary  of  the  upper  wing 
surface  begins  from  the  separation  point.  Furthermore,  an 
initial  direction  is  given  (such  as  the  median  value  between  the 
incoming  flow  direction  and  the  slope  of  the  wing  at  that  point) . 
Similarly,  the  lower  wing  surface  boundary  is  assumed  to  start 
from  the  trailing  edge.  An  initial  direction  is  also  given. 
Furthermore,  based  on  the  width  of  separation  at  the  trailing 
edge  (also  known  as  the  wake  altitude)  and  the  empirical  aspect 
ratio  of  the  wake; 


wake  length  L 

WF  =  wake  altitude  =  H  (See  Figure  3) 


it  is  possible  to  determine  the  initial  wake  length.  Moreover, 
both  upper  and  lower  boundaries  are  assumed  to  be  quadratic 
curves.  Thus,  the  initial  wake  boundary  can  be  determined.  The 
aspect  ratio  of  the  wake,  WF,  is  obtained  from  the  empirical 
curve  shown  in  Figure  4^1. 
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V.  Calculation  of  Boundary  Layer 


The  iteration  technique  only  considers  the  displacement 
thickness  effect  of  the  boundary  layer.  It  is  not  required  to 
study  the  structure  inside  the  boundary  layer.  Because  the 
integral  equation  method  is  relatively  simple  and  consumes  much 
less  time,  it  was  chosen  to  calculate  the  boundary  layer.  The 
calculation  of  the  boundary  layer  had  been  comprehensively 
introduced  in  reference  [1].  The  required  formula  and  method 
were  explained  in  detail.  The  following  is  a  brief  introduction 
of  this  method. 

(1)  Laminar  Boundary  Layer 

Due  to  the  fact  that  we  are  considering  incompressible 
flows,  the  Thwaite  method  was  chosen  to  calculate  the  laminar 
boundary  layer  on  the  wing. 

(2)  Estimation  of  Turning  Point 

The  Granville  method  was  adopted  to  determine  the  position 
of  the  turning  point.  Based  on  the  local  profile  factor  H  in 
laminar  boundary  layer  calculation,  we  used  the  criterion  curve 
derived  from  stability  theory  to  determine  whether  a  point  is 
unstable.  If  it  is  already  unstable,  it  is  necessary  to  use  the 
turning  criterion  curve  to  determine  whether  it  is  a  turning 
point.  Otherwise,  the  calculation  continues  based  on  laminar 
flow  condition. 

(3)  Turbulent  Boundary  Layer 

The  Nash  method  was  chosen  in  this  calculation.  This  is  a 
better  integral  equation  method  recommended  by  the 


two-dimensional  incompressible  turbulent  boundary  layer  meeting 
held  in  Stanford  in  1968.  This  method  considers  the  "upstream 
history  effect"  of  the  turbulent  boundary  layer. 

In  this  calculation,  we  assumed  that  separation  began  from 
H=»l .  8 . 

VI.  Results  and  Discussion 

In  reference  [1] ,  we  had  already  discussed  some  general 
problems,  such  as  the  effect  of  surface  element  division,  surface 
element  number  and  boundary  layer  calculation,  in  iterations.  In 
this  paper*,  the  simulation  of  separation  wake  at  large  attack 
angles  is  discussed.  In  addition,  two  classical  wings  are 
calculated.  The  computational  results  are  compared  to  the 
experimental  data. 

(1)  There  are  many  criteria  capable  of  determining  whether 

the  separation  wake  boundary  adjustment  in  internal  iteration  is 
convergent.  In  this  work,  it  is  determined  by  the  difference  of 
the  boundary  after  two  iterations.  It  was  found  that  the 
boundary  variation  becomes  very  small  after  two  adjustments.  The 
maximum  deviation  is  less  than  1%  chord  length.  Furthermore,  the 
maximum  is  located  downstream  far  away  from  the  wing.  Hence,  its 
effect  on  the  wing  surface  is  relatively  small.  Based  on  the 
fact  that  the  computed  results  of  both  wings  are  in  good 
agreement  with  the  experimental  data,  this  type  of  error  is  / 18 

certainly  allowable. 

(2)  The  rate  of  convergence  for  the  external  iteration  is 


also  very  fast.  At  large  attack  angles,  if  we  begin  to  calculate 
from  the  potential  flow  solution,  the  convergence  requirement  can 
generally  be  met  in  5~6  iterations;  i.e.,  after  two  consecutive 
iterations  the  change  in  c^  <  1%.  If  we  calculate  the  attack 
angle  continuously,  i.e.,  to  use  the  convergent  solution  of  the 
preceding  attack  angle  to  begin  to  solve  the  next  one,  a 
convergent  solution  can  usually  be  obtained  in  2~3  iterations. 

(3)  In  air  calculation,  it  was  found  that  the  initial 
direction  of  the  separation  wake  boundary  still  has  some  effect. 
At  small  attack  angles,  it  is  feasible  to  assume  that  it  is 
between  the  incoming  flow  direction  and  the  wing  slope  direction 
at  the  point  of  separation.  Nevertheless,  with  increasing  attack 
angle,  it  should  gradually  move  to  the  incoming  flow  direction  to 
make  the  computational  results  close  to  experimental  data. 

(4)  Figure  5  and  Figure  6  show  comparisons  of  experimental 
data  to  calculated  results  of  the  pressure  distribution  on  the 
Model  GA(W) -1  wing  when  the  attack  angles  are  16.04°  and  18.25°, 
respectively.  The  experimental  data  were  obtained  in  reference 
[8].  From  the  figures  we  can  see  that  they  agree  very  well  with 
the  exception  of  the  peak  values  at  the  leading  edge.  The 
position  of  separation  points  are  also  fairly  close.  Based  on 
comparisons,  the  experimental  data  to  calculated  results  when 
a»*i9.06*  and  20.05°  do  not  agree  as  well  as  at  other  attack 
angles.  In  this  case,  the  attack  angle  has  already  approached 
Clmax*  When  the  angle  of  attack  is  greater  than  20°,  the  forward 
movement  of  the  separation  point  on  the  upper  wing  surface  is 
very  large  and  sudden.  Therefore,  our  simple  model  cannot 
reflect  such  violent  changes  of  the  flow  pattern. 


K 


Figures  7  and  8  show  comparisons  of  the  calculated  and 


experimental  pressure  distributions  cf  the  NACA  4412  airfoil. 

The  experimental  data  were  obtained  from  reference  [9] .  Based  on 
these  figures,  they  agree  very  well.  The  experimental  data, 
however,  do  not  show  a  clear  position  of  the  separation  point  and 
the  separation  region.  This  might  be  due  to  measurement  errors 
because  based  on  force  measurements  separation  should  exist  on 
the  wing. 

/ 1 9 


1.  experimental 

2.  calculated 
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Figure  8 


1.  experimental 

2.  calculated 


(5)  Figures  9  and  10  show  the  experimental  data  and 
calculated  results  of  c^~a  and  cni~o  for  two  wing  models, 
respectively.  We  can  see  that  they  agree  well  until  Clmax.  This 
indicates  that  although  this  modeling  which  takes  separation 
effect  into  account  is  simple  and  coarse,  it  is  able  to  reflect 
the  reality.  However,  the  agreement  is  not  as  good  beyond  Clmax, 


which  indicates  that  this  model  is  limited. 


.-umm 


This  paper  presents  an  iterative  solution  for  calculating 
the  aerodynamic  force  on  a  wing  with  separation.  An  "equivalent 
body"  formed  by  streamline  displacement  is  used  to  replace  the 
original  geometric  body.  The  potential  flow  theory  is  used  to 
find  solution  for  this  equivalent  body.  It  is  capable  of 
obtaining  the  aerodynamic  characteristics  of  the  entire  wing 
until  near  Clmax.  The  boundary  of  the  separation  wake,  the  /20 

displacement  thickness  of  the  boundary  layer  and  the  position  of 
the  separation  point  are  iterated  as  a  part  of  the  whole 
solution.  Calculations  were  done  with  two  wing  models.  The 
results  show  that: 

(1)  It  is  feasible  to  use  this  streamline  displacement 
equivalent  body  concept  as  an  engineering  model  for  separation 
wake.  In  the  attack  angle  range  up  to  cl^^,  calculated  results 
agree  well  with  experimental  data. 

(2)  The  rate  of  convergence  of  this  computation  method  is 
very  fast,  especially  when  calculating  continuous  attack  angles. 

(3)  When  the  attack  angle  exceeds  clmax,  this  simple  model 
cannot  give  us  good  results.  It  remains  to  be  further  improved. 
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Numerical  Simulation  of  Vortex  Breakdown 
Shi  Xungang 
(Beijing  University) 

Abstract 

In  this  work,  the  breakdown  of  an  isolated  axisymmetric 
vortex  embedded  in  a  uniform  flow  is  examined  by  numerically 
integrating  the  complete  Navier-Stokes  equations  for  an  unsteady 
axisymmetric  flow.  The  results  show  that:  when  the  vortex  is 
small,  the  vortex  is  stable.  The  solution  approaches  a  steady 
flow.  A  quasi-cylindriCal  approximation  is  a  good  approximation 
for  this  type  of  flow.  When  the  vortex  strength ‘is  sufficiently 
high,  the  solution  is  unsteady.  A  re-circulation  zone  is  formed 
near  the  axis.  Its  shape  and  internal  structure  resemble  those 
of  the  multi-cell  breakdown  bubbles  observed  by  Faler  and 
Leibovich  (1978) .  With  appropriate  combinations  of  flow 
parameters,  the  flow  appears  to  be  quasi-per iodic  after  some 
time.  Parallel  calculations  made  with  the  quasi-cylindr ic 
approximation  show  that  as  far  as  predicting  the  breakdown  of  the 
vortex  is  concerned,  both  methods  are  in  good  agreement.  They 
both  show  that  in  the  relatively  low  Reynolds  number  range 
covered  by  this  work,  vortex  breakdown  is  not  closely  related  to 
Reynolds  number.  Furthermore,  it  is  independent  of  the  critical 
classification  of  the  upstream  flow. 


I.  Introduction 


Since  Peckham  and  Atkinson  (1957)  observed  the  breakdown  of 
leading  edge  vortex  of  a  sweepback  wing  at  a  large  attack  angle 
for  the  first  time,  various  theories  have  been  presented  by  a 
number  of  authors  to  explain  the  breakdown  of  the  vortex;  such  as 
the  critical  flow  theory  by  Squire  (1960)  and  Bossel  (1967, 

1969),  the  finite  transition  theory  by  Benjamin  (1962,  1967),  the 
fluid  dynamic  instability  theory  by  Ludwieg  (1962)  and  the  quasi- 
cylindrical  approximation  theory  by  Gartshore  (1962,  1963),  Hall 
(1965,  1966,  1967)  and  Mager  (1972).  Up  until  now,  however,  not 
a  single  theory  has  been  widely  accepted.  There  is  considerable 
confusion  among  various  theories  as  well  as  in  comparison  between 
theoretical  and  experimental  results. 

Based  either  on  the  critical  flow  theory  or  the  finite 
transition  theory,  the  f low  upstream  of  the  breakdown  must  be  super 
critical.  Hall  (1967,  1972)  and  Ludwieg  (1970)  also  pointed  out 
that  the  vortex  breakdown  process  described  by  quasi-cylindr ical 
approximation  is  a  process  approaching  the  critical  state  from  a 
super  critical  state.  Thus,  it  seems  to  be  a  relatively  popular 
viewpoint,  i.e.,  vortex  breakdown  must  begin  from  a  super 
critical  upstream  flow.  Nevertheless,  Lavan,  Nielsen,  and  Fejer 
(1969) ,  Kopecky  and  Torrance  (1973)  ,  and  Grabowski  and  Berger 
(1976)  numerically  integrated  the  Navier-Stokes  equations  for  a 
steady,  axisymmetric  flow  and  obtained  the  re-circulation  zone 
consisting  of  axisymmetric  breakdown  bubbles  similar  to  those 
observed  experimentally,  regardless  of  whether  the  upstream  flow  is 


super  critical  or  subcritical.  This  contradiction  still  remains 
to  be  classified. 


On  the  other  hand,  their  numerical  results  failed  to  show 
the  double-cell  structure  inside  the  breakdown  bubble  as  observed 
by  Paler  and  Leibovich  (1978)  by  laser  flowmeter  (See  Figure  1) . 
Furthermore,  when 
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the  rotational  velocity  is  sufficiently  high,  their  calculations  /23 
could  not  even  result  in  convergent  solutions.  This  might  be 
attributed,  to  the  fact  that  the  assumptions  they  made  are  too 
rigorous.  Leibovich  (1978)  believed  that  numerical  results  could 
not  describe  such  double  cell  structure  if  the  experimentally 
observed  periodicity  and  asymmetry  of  the  flow  inside  the 
breakdown  bubbles  were  not  considered. 

In  this  work,  their  opinions  are  noted.  Through  numerical 
integration  of  the  complete  N-S  equations  for  an  unsteady 
axisymmetr ical  flow,  we  wish  to  be  able  to  more  closely  describe 
the  breakdown  of  an  "axisymmetr ic"  vortex.  Based  on  economic 
consideration  and  computer  limitation,  the  axisymmetr ic 
assumption  is  reluctantly  kept.  In  order  to  facilitate  the 
analysis  and  comparison,  parallel  calculations  based  on  quasi- 
cylindrical  approximation  are  also  made.  On  this  basis,  the 
contradiction  between  theoretical  and  numerical  results  mentioned 
above  is  discussed. 
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II.  Mathematical  Model 


In  this  paper,  an  isolated  axisymmetric  vortex  embedded  in 
an  unbounded  uniform  flow  is  studied.  Let  us  choose  a 
cylindrical  coordinate.  The  x-axis  coincides  with  the  axis  of 
symmetry.  Let  us  assume  that  the  axial  and  circumferential 
velocity  distributions  on  an  "inlet  cross-section"  are  given  by 
the  two  following  equations: 

x  =0t  {  “=!+«/('> 

1  w  —T ,  g(r) 


where  the  .velocity  components  and  coordinates  are  rendered 
dimensionless  by  using  the  velocity  of  the  uniform  flow  at 
infinity  u«>  and  vortex  core  radius  R.  Figure  2  schematically 
shows  their  profiles.  The  initial  axial  velocity  on  the  axis, 
Uq=1+«.  and  the  velocity  at  infinity  or  the  vortex  strength  Tg 
are  the  two  flow  parameters. 


Figure  1.  Flow  Field  Structure  Inside  an  Axisymmetric  Breakdown 
Bubble  (Leibovich  1978) 
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Figure  2.  Initial  Axial  and  Circumferential  Velocity 
Distributions 


If  we  introduce  a  local  moment,  r=  r w 
the  circumferential  vortex  component  dv/ax  -  3u/ar 
and  the  flow  function  V, 

then  the  dimensionless  N-S  equation  for  an  unsteady,  axisymmetric 
flow  of  a  viscous,  incompressible  fluid  can  be  re-written  as  a 
series  of  equations  of  f-Q-W. 
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where  Re-umRlv  This  is  a  series  of  parabolic  equations.  In 
terms  of  spatial  coordinates  alone,  it  is  also  elliptical.  In 
order  to  obtain  a  unique  solution,  it  is  necessary  to  have  the 
appropriate  initial  conditions  and  the  boundary  conditions  along 
the  entire  boundary  of  integration. 

The  integration  region  is  given  as  follows: 

D={q^.x<,L,  0<r<  +  00} 

The  outlet  end  of  the  boundary  is  chosen  sufficiently  for 
downstream  at  x=L,  where  L»l.  For  example,  we  choose  L=20 

The  conditions  to  have  a  unique  solution  are  as  follows: 
Initial  condition: 

*-0,  r-r,rg(r),  Q  =  -af\r),  V- -«f' r/<r)</r  + 

The  boundary  conditions  in  D  are: 
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The  so-called  quasi-cylindrical  approximation  is  to  assume 


V  «u,  w 


Furthermore,  the  flow  is  steady.  In  this  approximation,  the 
above  equations  can  be  simplified  as: 
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This  is  also  a  series  of  parabolic  equations.  After  the  initial 
conditions  at  the  initial  cross-section  x=0  and  the  boundary 
conditions  at  the  axis,  r=0,  and  at  the  outer  edge  r  -+>+<»,  are 
given,  it  is  always  possible  to  obtain  a  numerical  solution  along 
the  x  direction.  The  initial  conditions  are: 

x-0 :  r=r0rg(r),  D=-af'(r) 

The  boundary  conditions  are: 

r=0  f=  0  0  =  0 


r-*.+«a 


r— r  Q-+o 
1  l0 


However,  the  quasi-cylindrical  approximation  is  no  longer 
valid  near  the  vortex  breakdown  point.  The  differential 
equations  also  become  unsteady.  Numerical  calculations  are  no 
longer  convergent  from  this  point  onward.  Hence,  the  appearance 
of  large  axial  gradient  and  computational  divergence  can  be 
considered  as  a  sign  for  vortex  breakdown. 

In  order  to  make  the  infinite  integration  area  in  the  r 
direction  finite,  and  to  ensure  that  the  numerical  solution  has 
sufficiently  high  resolution  in  areas  where  the  flow  varies 
vigorously,  two  independent  coordinate  transformations  are 
introduced  for  the  radial  and  axial  direction,  respectively: 
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In  the  quasi-cylindrical  approximation,  only  the  radial 
transformation  is  required.  The  above  differential  equations,  as 
well  as  the  initial  and  boundary  conditions,  must  be  transformed 
correspondingly. 

III.  Computational  Results  and  Discussion 

The  average  implicit  finite  difference  method  developed  by 
Crank-Nicolson  was  used  to  solve  the  simplified  equations  in 
quasi-cylindrical  approximation.  This  technique  is  simple  and 
efficient.  The  results  show  that  axial  velocity  varies  very 
slightly  along  the  axis  when  Tg  is  very  small,  e.g.,  r0=0.63.  If 
we  continue  to  proceed  downstream,  axial  velocity  decreases 
initially  and  then  rises  gradually;  approaching  the  incoming  flow 
velocity  U-+1.  The  vortex  is  stable.  As  Tg  increases,  axial 
velocity  decreases  more  rapidly.  When  Tq  drops  below  a  certain 
value,  axial  velocity  drops  abruptly  to  a  value  very  close  to  0 
in  front  of  a  specific  position.  The  calculation  is  no  longer 
convergent  afterward.  As  mentioned  above,  this  position  can  be 
considered  as  the  breakdown  point.  With  increasing  Tg,  the 
breakdown  point  position  continues  to  move  upstream.  An  increase 
in  Og,  however,  will  strengthen  the  stability  of  the  vortex.  The 
sumbols  in  Figure  3  represent  a  summary  of  calculated  results 
obtained  using  the  quasi-cylindrical  approximation.  Based  on  the 


plot,  with  the  exception  of  one  point  Uq=1.4  and  ^=0.8944,  the 
conclusions  on  whether  the  vortex  breaks  down  or  not  agree  with 
those  obtained  by  calculation  with  Re=100  and  200.  The  dot  and 
line  curve  in  the  figure  represents  the  vortex  breakdown 
boundary. 


Figure  3.  Summary  of  Calculated  Results 

1.  stable 

2.  unstable 

3.  super  critical  region 

4.  subcritical  region 


An  alternate  direction  iteration  method  was  used  to  solve  /26 
the  complete  N-s  equations.  The  results  show  that  a  stable  flow 
predicted  by  the  quasi-cylindrical  approximation  reaches  a 
certain  steady  state  after  some  time.  The  flow  surface  appears 
to  be  very  smooth.  The  quasi-cylindrical  approximation  should 
obviously  be  valid  in  the  entire  flow  field.  A  comparison  of 
velocity  distribution  also  shows  that  the  flow  fields  obtained  by 
two  different  methods  are  in  good  agreement.  This  indicates  that 
the  quasi-cylindrical  approximation  is  appropriate  for  a  stable 
vortex.  In  addition,  the  computer  program  designed  for 
numerically  integrating  the  complete  N-S  equations  also  underwent 
a  thorough  test. 

When  predicting  vortex  breakdown  under  quasi-cylindrical 
approximation,  except  for  individual  edge  cases,  the  numeric 
solution  of  the  complete  N-S  equations  shows  that  the  flow  is 
unsteady.  Figure  4  is  an  example.  They  are  a  series  of  cross- 
sections  of  an  axisymmetric  flow  surface  and  a  meridian  plane. 

As  time  goes  by,  it  first  begins  to  bulge  near  the  axis  and  then 
develops  into  an  enclosed  re-circulation  zone.  This  re¬ 
circulation  zone  continues  to  grow  into  a  so-called  "double-cell" 
structure  (See  Figure  4*)  .  Its  shape  and  internal  structure  are 
very  similar  to  those  of  the  breakdown  bubble  obtained  by  Faler 
and  Leibovich  (1978)  with  a  laser  flowmeter  (See  Figure  1). 

Under  certain  appropriate  combinations  of  flow  parameters,  the 
flow  exhibits  some  quasi-periodicity  after  some  time.  A  new 
inner  cell  is  continuously  formed  periodically  at  the  tip  of  the 
breakdown  bubble.  It  gradually  intensifies  and  moves  backward 
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The  periodicity  of  the  unstable  vortex  flow  in  our  numerical 


results  is  qualitatively  in  agreement  with  the  experimental 
observation  made  by  Sarpkaya  (1971) ,  and  Faler  and  Leibovich 
(1978)  . 

If  we  assume  that  the  flow  becomes  unstable  with  respect  to  / 28 
any  non-ax i symmetric  disturbance  after  the  breakdown  of  the  "axi- 
symmetric"  bubble,  then  the  inner  cells  which  periodically  flow 
downstream  might  be  the  spiral  wake  following  the  breakdown  of  a 


K 


bubble  as  observed  in  the  experiment. 

In  agreement  with  the  conclusions  obtained  with  the  quasi- 
cylindrical  approximation,  when  the  complete  N-S  equations  are 
used  in  the  calculation,  as  far  as  whether  the  vortex  breaks  down 
or  not,  the  results  totally  coincide. with  those  calculated  based 
on  Re=100  and  200.  The  solid  line  in  Figure  3  represents  the 
vortex  breakdown  boundary  calculated  based  on  the  complete  N-S 
equations  for  an  unsteady  flow.  It  is  between  the  curve  obtained 
based  on  the  quasi-cylindrical  approximation  and  the  curve 
(dotted)  obtained  by  Grabowski  and  Berger  (1976)  based  on  the  N-S 
equation  for  a  steady  flow.  However,  they  are  very  close. 

Figure  3  also  shows  the  critical  curve  (double  dotted  line) 
separating  super  critical  upstream  flow  from  subcritical  upstream 
flow  as  obtained  by  Mager  (1972)  based  on  Benjamin's  checking 
equation  (1962) .  Based  on  the  figure,  just  as  Grabowski  and 
Berger  (1976)  pointed  out,  many  breakdown  solutions  were  obtained 
with  subcritical  upstream  conditions.  It  is  not  in  agreement 
with  Squire's  critical  flow  theory  and  Benjamin's  finite 
transition  theory.  Parallel  calculations  based  on  the 
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quasi-cylindrical  approximation  also  proved  that  these 
subcritical  upstream  flows  would  definitely  lead  to  vortex 
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breakdown.  The  parabolic  nature  of  the  quasi-cylindrical 
approximation  equations  eliminates  the  possibility  of  the 
propagation  of  any  disturbance  upstream.  Hencef  the  vortex 
breakdown  thus  calculated  could  not  be  explained  based  on  the 
upstream  propagation  of  disturbance.  Several  experiments 
conducted  in  pipes  also  indicate  that  when  the  flow  remains 
unchanged  and  the  rotation  is  strengthened,  as  the  subcritical 
nature  of  the  upstream  flow  is  strengthened  based  on  Mager's 
critical  curve,  only  the  breakdown  point  is  moved  up  and  the 
pattern  is  changed.  It  does  not  mean  that  the  breakdown 
disappears.  Therefore,  the  only  conclusion  we  could  reach  is 
that  vortex  breakdown  is  indeed  independent  of  the  subcritical 
nature  of  the  upstream  flow. 
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IV.  Conclusions 


In  this  paper  the  breakdown  of  an  isolated  axisymmetric 
vortex  embedded  in  a  uniform  flow  was  investigated  through  the 
numeric  integration  of  the  N-S  equations  at  relatively  low 
Reynolds  numbers. 

First,  the  N-S  equations  were  simplified  by  the  quasi- 
cylindric  approximation.  Its  solution  can  be  numerically 
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determined  step  by  step  along  the  x-direction.  The  results  show 
that  it  is  an  excellent  approximation  of  the  actual  flow  for  a 
stable  vortex.  Any  sudden  deceleration  of  the  axial  flow  and 
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divergence  of  the  calculation  can  be  considered  as  a  sign  for 
vortex  breakdown.  Furthermore,  vortex  breakdown  is  very 
sensitive  to  changes  in  vortex  strength. 

Then,  the  N-S  equations  for  an  unsteady  axisymmetr ic  flow 
are  numerically  integrated.  If  vortex  breakdown  does  not  take 
place,  the  solution  will  approach  a  steady  state.  Otherwise,  the 
solution  will  remain  unsteady.  A  re-circulation  zone  will  appear 
near  the  axis.  Its  shape  and  internal  structure  are  very  similar 
to  those  of  the  breakdown  bubble  observed  experimentally  by  Faler 
and  Leibovich  (1978) .  Under  appropriate  combinations  of  flow 
parameters,  the  flow  will  appear  periodic  after  some  time. 

The  results  of  both  methods  show  that  as  far  as  the 
relatively  low  Reynolds  numbers  covered  by  our  calculations  are 
concerned,  vortes  breakdown  is  not  very  closely  related  to  Reynolds 
number.  Moreover,  it  is  independent  of  the  critical 
classification  of  the  upstream  flow. 
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Environmental  Study  on  Turbulent  Boundary  Layer  Separation  / 30 

Wei  Qingding 
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Introduction 

Boundary  layer  separation,  as  well  as  effects  such  as 
stalling  and  choking,  are  very  important  problems.  This  is  not 
only  because  great  difficulties  and  losses  of  aircraft  and 
various  fluid  equipment  may  be  the  result  of  stalling 
and  choking,  but  also  because  such  equipment  is  often  in  its 
optimum  operating  condition  prior  to  stalling  and  choking. 

Therefore,  it  is  very  important  to  clearly  grasp  the  mechanism  of 
flow  separation  and  to  effectively  predict  and  control 
separation.  It  is  a  very  practical  topic.  Over  the  years,  many 
researchers  dedicated  their  work  in  this  area.  Considerable 
progress  has  been  made  in  theory  and  application  as  well.  The 
successful  application  of  boundary  layer  control  on  aircraft  is 
an  example.  Nevertheless,  we  cannot  say  that  we  clearly 
understand  boundary  layer  separation,  especially  the  separation 
of  a  turbulent  boundary  layer.  Many  studies  just  begin  to  become 
more  in-depth.  Due  to  the  unsteadiness  of  turbulence  and 
difficulties  in  dealing  with  three-dimensional  and  random 
problems,  there  are  special  difficulties  in  treating  theoretical 
and  experimental  aspects  of  turbulent  boundary  layer  separation. 

In  terms  of  experimental  research,  there  was  almost  no  effective 
means  to  study  the  flow  structure  of  turbulent  boundary  layer 


separation  before  the  70' s.  In  the  early  50' s,  although 
Schubauer,  G.B.  and  Klebanoff,  P.S.t1!  did  some  fine  work  using 
available  measuring  techniques,  the  understanding  gained  on  the 
boundary  layer  was  limited  to  the  classical  model.  Not  only  the 
changing  nature  of  turbulent  boundary  layer  separation  was  not 
involved,  but  also  the  flow  structure  in  the  separation  zone  was 
touched.  As  flow  visualization  techniques,  hot  wire  techniques 
and  laser  technique  made  progress,  it  was  clear  in  the  60' s  that  a 
turbulent  boundary  layer  is  very  different  from  a  laminar  boundary 
layer.  Laminar  flow  separation  is  always  along  a  clear  line-a 
two  or  three  dimensional  separation  line.  Turbulent  flow 
separation,  however,  is  a  transient  process,  from  weak  to  strong. 
Sandborn,  V.A.  and  Kline,  S.J.  (1961)  ^  called  the  former  steady 
or  defined  separation  and  the  latter  the  changing  or  intermittent 
separation.  Nevertheless,  they  did  not  quantitatively  describe 
the  transient  nature  of  turbulent  flow  separation.  After  the 
70' s,  many  researchers  used  techniques  such  as  flow 
visualization,  hot  line  and  laser  to  do  a  great  deal  of  delicate 
studies  on  the  turbulent  boundary  layer  separation  zone  to  bring 
the  experimental  work  on  turbulent  boundary  layer  into  a  new 
quantitative  stage.  Simpson ^ 3—7 ]  et  an(j  Kline^®-*^  et  al 
have  systematically  worked  in  this  area.  Simpson,  in  particular, 
used  laser,  hot  wire,  and  hot  film  techniques  to  systematically 
measure  various  flow  parameters  in  the  turbulent  boundary  layer 
separation  zone.  The  author,  under  the  guidance  of  Professor 
Hiroshi  Sato^“^  ,  began  to  conduct  a  delicate  investigation  on 
the  turbulent  boundary  layer  separation  zone  by  using  hot  wire 


techniques  such  as  a  compound  probe  for  wind  speed  and  direction 
as  well  as  by  using  a  smoke  flow  visualization  technique.  In 
particular,  the  effect  of  velocity  distribution  of  the  main  flow 
on  turbulent  boundary  layer  separation  was  investigated  based  on 
dynamic  characteristics.  In  addition,  the  effect  of  large  scale 
motion  upstream  on  the  intermittent  separation  of  turbulent 
boundary  layer  was  also  studied.  This  paper  briefly  describes 
our  work.  In  order  to  facilitate  the  comparison,  the  study  done 
by  Simpson  et  al  is  also  appropriately  introduced. 

Manuscript  received  on  June  20,  1984 

II.  Primary  Methods  in  Studying  Turbulent  Boundary  Layer  / 3 1 

Separation 

1.  Experimental  Apparatus 

The  experimental  apparatus  should  be  able  to  fully  develop  a 
turbulent  boundary  layer  and  cause  its  separation  at  a  negative 
pressure  gradient.  Our  work  was  conducted  in  the  lmximxsm 
experimental  section  of  a  low  velocity  wind  tunnel  as  shown  in 
Figure  1.  Due  to  the  tripping  wire,  the  flow  creates  a  turbulent 
boundary  layer  over  the  plate  which  separates  on  the  parabolic 
surface  because  of  negative  pressure,  upstream  the  plate,  there 
are  devices  to  create  velocity  and  temperature  cross-sections. 

The  apparatus  used  by  Simpson  was  a  wind  tunnel  with  a  curved  top 
plate  such  as  the  one  shown  in  Figure  2.  The  curved  wall  was 
designed  to  create  a  pressure  distribution  similar  to  that  on  the 
wing  cross-section.  The  turbulent  boundary  layer  was  formed  by 


the  tripping  effect  of  the  dull  leading  edge  of  the  plywood. 
These  two  specific  devices  and  other  published  experimental 
apparatus  are  limited  to  the  two-dimensional  case. 
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Figure  1.  Longitudinal  Cross-section  of  the  imMmxsm  Low 
Velocity  Wind  Tunnel  and  Its  Coordinate  System 

1.  tripping  wire 


Figure  2.  Schematic  Diagram  of  the  Experimental  Section  Used  By 
Simpson  to  Study  Turbulent  Boundary  Layer.  Each 
major  division  is  10  inches. 
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2.  Major  Testing  Equipment 

Due  to  the  complexity  of  the  turbulent  separation  flow,  the 
measuring  apparatus  must  have  a  specific  frequency  response 
(ordinary  frequency  response  at  10  KHz)  and  directivity  (at  least 
capable  of  recognizing  forward  reverse  flows) .  Its  interference 
to  the  flow  field  must  be  held  at  a  minimum.  Any  conventional 
hot  wire  tachometer  even  equipped  with  an  X  probe  cannot  meet  the 
directivity  requirement.  Simpson  et  al  used  laser  to  aid 
compound  hot  wire  flow  meter  and  a  three  wire  probe  as  shown  in 
Figure  3  to  make  the  measurement.  This  method  is  simpler  than 
using  a  laser  in  probe  installation,  spatial  measurement  and  data 
processing.  The  cost  is  also  lower.  However,  it  cannot  measure 
as  close  to  the  wall  as  the  laser.  Laser  can  measure  as  close  as 
0.0015,  while  a  three  wire  probe  can  only  measure  to  0.025. 

No  matter  whether  a  laser  or  a  hot  wire  is  used,  a  computer  / 32 
is  needed  to  process  the  data  rapidly.  In  our  experiment,  we 
used  an  on-line  real  time  data  processing  system. 


Figure  3.  Three-wire  Wind  Speed  Wind  Ditection  Probe 

1.  mm 

2.  mm 

3 .  mm 

Flow  visualization  techniques  are  important  in  studying  the 
separation  of  a  turbulent  boundary  layer.  Due  to  the  fact  that  a 
turbulent  flow  field  is  a  random  three-dimensional  unsteady  flow 
field,  a  laser  or  hot  wire  flow  meter  which  measures  a  point  can 
only  provide  some  instantaneous  or  statistical  information  at 
that  point.  A  flow  visualization  technique,  however,  can  be  used 
to  obtain  the  flow  pattern  of  the  entire  field.  Accurate  flow 
visualization  results  can  also  provide  quantitative  information. 
In  our  experiment,  we  used  a  smoke  method.  In  addition,  we  used 
a  multi-channel  short  exposure  photography  technique  and  a 
controllable  condition  method  to  improve  the  objectivity  and 
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controlability  (to  select  characteristic  flow  patterns)  of  flow 
visualization. 

III.  Major  Experimental  Results 

1.  Statistical  Quantity  Describing  Intermittent  Separation 
of  Turbulent  Boundary  Layer. 

The  separation  of  a  turbulent  boundary  layer  is  intermittent. 
In  the  separation  area,  separation  and  attachment  take  place 
alternately.  This  can  be  reflected  by  the  intermittent  reverse 
flow  in  the  flow  extremely  close  to  the  wall  surface.  The 
probability  of  back  flow  represents  the  degree  of  separation. 

Let  us  define  a  statistical  average-back  flow  rate  Rn, 

where, 

r,-2 

r.= 2  a/..  o<  t  <r 

•  _ 

Atpi  and  Atni  are  the  time  intervals  over  which  forward  and 
backward  flow  occur  at  the  point  of  measurement,  respectively. 
They  are  determined  by  the  forward  and  backward  flow  signals 
obtained  by  the  three-wire  wind  speed  wind  direction  probe  (See 
Figure  4) . 

Based  on  the  frequency  histogram  measured  with  laser, 
Simpson  defined  a  similar  quantity  Vp.  It  represents  the 


probability  for  forward  flow.  When  there  is  not  backward  flow, 


Figure  5  shows  the  backward  flow  rate  distribution  measured 
when  the  incoming  flow  is  a  3m/s  uniform  flow.  Figure  6  shows 
Simpson's  results  based  on  an  18m/s  uniform  flow.  These  results 
are  similar.  There  is  no  back  flow  upstream,  i.e.,  Rn=0,  Vp=l. 

At  certain  spots  near  the  wall  (x«20cm  in  Figure  5  and  x«122.0  in 
Figure  6} ,  intermittent  backward  flows  of  low  probabilities  begin 
to  appear.  Further  downstream,  the  probability  of  backward  flow 
increases  (Rnf,  Kp|) •  The  backward  flow  area  in  the  vertical 
wall  direction  becomes  thicker.  In  the  deep  separation  region 
(around  x=50cm  in  Figure  5  and  x>150  in.  in  Figure  6) ,  the 
probability  for  backward  flow  to  occur  is  greater  than  0.90. 
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Figure  4.  Forward  Flow  Signal  Sp,  Backward  Flow  Signal  Sn, 
and  0-1  Square  Wave  Defined  by  Respective  Thresh 


Values. 


by  Respective  Threshold 


1.  threshold  value 

2.  threshold  value 


TOW 


Figure  5.  Distribution  of  Backward  Flow  Rate  Rn  at  Various 

Cross-sections  in  the  Separation  Area  when  the  Main 
Flow  Wind  Speed  is  3m/s. 


Figure  6.  Forward  Flow  Rate  V_  Distribution  in  Turbulent 

Boundary  Layer  Separation  Region  Measured  by  Simpson. 

(a)  Distribution  of  }>  in  x-direction  along  wall 
surface  at  0 . 25iturr  height 

(b)  Distribution  of  Vp  in  five  positions 

1.  position  (in.) 

2.  Flow  Visualization  of  Separation  Area 

Figure  7  shows  the  photograph  of  the  flow  in  the  separation 
region  and  its  upstream  area  as  shown  by  the  smoke  line  method. 
The  smoke  lines  are  parallel  to  the  wall  and  perpendicular  to  the 
flow.  The  smoke  track  shows  the  instantaneous  wind  speed 
distibution  6mm  (0.075§)  away  from  the  wall.  We  can  see  that 
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back  flow  begins  to  occur  at  x=12cm.  There  is  not  backward  flow 
upstream  from  this  position.  Downstream  from  this  position,  the 
back  flow  region  begins  to  expand  gradually.  This  phenomenon  is 
in  agreement  with  that  measured  by  a  wind  speed  wind  direction 
meter.  Figure  8  shows  the  instantaneous  backward  flow  region 
shown  in  Figure  7.  We  can  see  the  transitional,  intermittent  and 
three-dimensional  nature  of  turbulent  boundary  layer  separation. 

3.  Average  Velocity  Distribution  in  Separation  Region 
In  the  boundary  layer  separation  region,  the  direction  and 
magnitude  of  the  instantaneous  velocity  at  every  point  vary 
continuously.  The  average  speed  we  are  measuring  is  the 
algebraic  average  of  the  forward  and  backward  flow. 

2  U,(x,y) 


where  N  is  the  number  of  sampling  points  and  U^(x,y)  is  the 
component  of  the  instantaneous  wind  velocity  on  the  x-y  plane  as 
measured  by  a  wind  speed  wind  direction  meter  which  is  positive 
for  forward  flow  and  negative  for  backward  flow. 

Figure  9  shows  the  average  wind  velocity  distribution.  As 
we  can  see  that  the  backward  flow  rate  is  zero  at  curvature 
varying  points  on  an  average  velocity  cross-section.  At  points 
where  the  average  velocity  is  zero,  the  backward  flow  rate  is 
approximately  0.50.  Figure  10  shows  similar  results  obtained  by 
Simpson. 
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Figure  7.  Flow  Visuali2ation  of  Turbulent  Boundary  Layer 

Separation  Before  and  After  Separation  (Smoke  Method* 
Smoke  yw=0.075 5  f rom  Wall) 


1.  curved  wall 

2.  (cm) 

3.  20cm 


T3 
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Figure  8.  Instantaneous  Separation  Region  in  Turbulent  Boundary 
Layer 

.  1.  leading  edge  of  separation  area 

2.  (cm) 

3.  backward  flow  area 
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Figure  9. 


Distribution  of  Average  Velocity,  Turbulence  and  Back 
Flow  Rate  in  Turbulent  Boundary  Layer  Separation 
Region  (U^=3m/s) 
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4.  Pulse  Velocity  Distribution  in  Separation  Region 
The  curve  shown  in  Figure  9  shows  the  typical  turbulence 


distribution  measured  in  separation  area.  As  backward  flow 
begins  to  appear,  the  corresponding  velocity  pulse  is  the 
maximum.  Near  the  border  of  the  backward  flow  region  (Rn=0.01), 
velocity  pulse  is  also  the  highest.  Thus,  the  dimensionless 
pulse  velocity  is  determined  by  the  following  equation: 


u '  (x.y)/l’, 


U  (x.ij) ) 


l 


where  is  the  incoming  flow  velocity. 
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Figure  10.  Mean  Wind  Speed  Distribution  in  Turbulent  Boundary 
Layer  Separation  Region  Measured  by  Simpson  Using 
Pitot  Tube  (+) ,  Hot  Film  (0)  and  Laser  (A,n) 
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Figure  12.  Wind  Direction  on  xz  Plan^in  Turbulent  Boundary 

Layer  Separation  Zone  (x=55,y  corresponds  to  a.  1.0; 
b.  3.0;  c.  5.0;  d.  7.0;  e.  11.0;  and  £  15.0  (cm), 
respectively) 


5.  Wind  Direction  Probability  Distribution  in  Separation  / 36 

Zone 

We  used  5-wire  and  8-wire  probes  to  measure  the  wind 
direction  probability  distribution  in  the  separation  zone. 

Figures  11  and  12  show  the  wind  direction  distributions  on  x-y 
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and  x-z  planes,  respectively.  We  can  see  that  in  the  deeply 
separated  region,  the  wind  direction  distribution  range  covers 
the  entire  sphere.  Around  Rn=0.5,  the  probability  in  any 
direction  is  very  close. 

6.  Effect  of  Main  Flow  Velocity  Distribution  on  Turbulent 
Boundary  Layer  Separation 

Figure  13  compares  the  intermittent  separation  zones  which 
correspond  to  five  main  flow  velocity  distributions  (uniform 
flows  at  lm/s,  3m/s  and  5m/s  as  well  as  uniform  positive  and 
negative  flows  at  3m/s  at  the  boundary  layer) .  It  is  possible  to 
see  that  the  separation  zone  corresponding  to  the  negative  shear 
main  flow  is  the  smallest.  The  separation  zone  corresponding  to 
the  main  flow  with  a  positive  shear  flow  is  the  largest.  For 
uniform  flows,  the  separation  zone  becomes  smaller  with 
increasing  Re. 


a  13 


Figure  13.  Turbulent  Boundary  Flow  Separation  Zones  Correspond¬ 
ing  to  Five  Main  Flow  Velocity  Distributions 
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1.  1 (m/sec)  2.  3 (m/sec) 

3.  5 (m/sec)  4.  uniform  main  flow 

5.  positive  shear  main  flow 

6.  negative  shear  main  flow 


Figure  14.  Apparatus  Used  to  Measure  the  Correlation  Between 
Upstream  Large  Scale  Movement  and  Intermittent 
Separation  in  the  Front  Section  of  the  Separation 
Zone 

1.  single  hot  wire  probe 

2.  triple  hot  wire  wind  speed  and  direction  probe 

3.  separation  zone 


7.  Effect  of  Large  Scale  Movement  Upstream  on  Intermittent 
Separation  of  Turbulent  Boundary  Layer 

Figure  14  is  a  schematic  diagram  of  the  apparatus  used  to 


measure  the  correlation  between  any  large  scale  movement  upstream 
from  the  separation  zone  and  the  intermittent  separation  in  the 


front  section  of  the  separation  zone.  A  single  hot  wire  probe 
was  installed  upstream.  The  velocity  signal  measured  is  U(t).  A 
wind  speed  and  direction  probe  was  installed  downstream  in  the 
separation  zone.  The  backward  flow  signal  to  be  detected  is 
Sn(t) .  A  conditional  sampling  method  was  used.  The  time  t^  at 
which  a  backward  flow  appears  might  be  used  as  the  base  to 
measure  the  velocity  signal  in  the  time  interval  [t^-^r^f  t^+^Tj] 
(See  Figure  15) .  Figure  16  shows  an  example  of  the  result  of  the 
measured  conditional  average  speed.  Based  on  the  conditional 
sampling  average  velocity,  it  is  possible  to  obtain  the  velocity 
cross-section  changes  in  a  specific  time  interval.  Figure  17 
shows  the  .variation  of  the  conditional  average  velocity  on  seven 
cross-sections  along  the  straight  and  curved  walls  at  seven  time 
intervals  within  t*-ll/30~0  sec.  The  solid  lines  represent  the 
conditional  average  velocity  and  the  dotted  lines  are  the  usual 
average  velocity.  The  velocity  in  slanted  line  covered  areas  are 
lower  than  average.  They  are  called  "low  velocity  flow  blocks". 


Figure  15.  Sampling  Conditional  Signal  Sn(t)  and  Conditional 
Sampling  Signal  U(t) 

Areas  covered  by  crossed  slanted  lines  are  higher  than  average 
velocity  areas  -  called  "high  velocity  flow  blocks".  Based  on  / 37 
the  figure,  as  time  progresses,  a  low  velocity  flow  block  flows 
downstream  following  a  high  velocity  flow  block.  When  the  low 
velocity  block  flows  by  the  intermittent  separation  zone, 
backward  flow  occurs.  We  can  thus  conclude  that  any  large  scale 
movement  upstream  can  affect  intermittent  separation.  Low 
velocity  flow  block  is  one  of  the  causes  for  inducing 
intermittent  separation. 


Figure  16.  An  Example  of  Conditional  Average  Wind  Velocity 

1.  0(t)  (0.02  m/sec) 

2.  probe  position  (x,y  cm) 

3.  t(sec) 
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Figure  17.  Experimental  Results  on  the  Correlation  Between  the 

Low  Velocity  Flow  Block  Upstream  and  the  Intermittent 
Variation  of  the  Leading  Edge  of  the  Separation  Zone 


IV.  Major  Conclusions  and  Discussion 


The  separation  of  a  turbulent  boundary  layer  is  intermittent, 
transitional  and  three-dimensional.  It  is  possible  to  use 
forward  and  reverse  flow  rates  to  describe  the  extent  of 
separation.  The  distributions  of  certain  statistical  averages 
such  as  average  velocity  and  pulse  velocity  is  related  to  the 
distribution  of  backward  flow  rate.  The  velocity  distribution  of 
the  main  flow  will  determine  the  ease  of  separation.  Large  scale 
upstream  motion  will  affect  the  intermittent  separation  of  tne 
turbulent  boundary  layer. 

The  primary  objective  to  study  turbulent  boundary  layer 
separation  is  to  predict  and  control  boundary  layer  separation. 
The  purpose  of  the  experiment  is  to  clarify  the  physical  nature 
of  the  flow  phenomenon,  and  to  provide  effective  theoretical 
models  or  empirical  or  semi-empirical  formulas.  Due  to  the 
complexity  of  the  problem,  these  tasks  are  not  yet  accomplished. 
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Turbulent  Boundary  Layer  Characteristics  In  and  Out  of 
Separation  Region  Measured  at  Wing-Plate  Junction  at  Low  and 

Subsonic  Speeds 
Xin  Dingding  and  Deng  Xueying 
(Beijing  Institute  of  Aeronautics  and  Astronautics) 

Abstract 

This  paper  presents  the  three-dimensional  turbulent  boundary 
layer  measurement  made  in  and  out  of  the  separation  region  at  the 
wing-plate  junction  at  low  and  subsonic  speeds.  In  the  low  speed 
experiment-,  the  distributions  of  the  three  average  velocity 
components  and  six  turbulent  stress  components  at  ten  stations 
along  two  lines  were  measured.  At  subsonic  speeds,  similar 
measurements  were  made  at  two  stations  in  and  out  of  the 
separation  region.  The  results  show  that  separation  induces  the 
formation  of  a  horseshoe  shaped  vortex,  which  is  an  important 
factor  affecting  the  behavior  of  the  three-dimensional  boundary 
layer  in  and  out  of  the  separation  region.  In  addition,  method 
to  process  experimental  data  is  also  discussed.  Experimental 
results  influenced  by  pressure  gradient  and  streamline  curvature 
were  also  analyzed. 


Symbols 


A,B,C,D,E,F  coefficients  in  equation  (6)  and  (7) 
c,h  coefficients  in  equation  (3) 


Ec  voltage  output  from  linearizer 

ec  pulse  of  voltage  output  from  linearizer 

K  hot  wire  calibration  index  in  equation  (2) 

M  local  Mach  number 

Moo  incoming  flow  Mach  number 

Am  deviation  in  Mach  number  distribution 

P  static  pressure  of  flow 

T  static  temperature  of  flow 

u  local  flow  velocity 

u.o  incoming  flow  velocity 

ue  velocity  at  outer  edge  of  boundary  layer 

ue££  •  cooling  rate  of  hot  wire  in  equation  (2)  and  (3) 

ux,uy,u2  velocity  components  along  (x,y,z)  axes 

us,uyrun  velocity  components  along  (s,y,n)  axes 

UT,UN,UB  velocity  components  along  (T,N,B)  axes 

(u'.u,).  =  (pV.O/P2 »  mass  weighted  shear  stress 

a,  angles  indicating  the  orientation  of  the  hot  wire 

(See  Figure  3) 
p  density  of  flow 

superscript 

(-)  time  average  or  mean  value 

(')  fluctuation 

subscript 

A,B  output  signals  from  hot  wires  A  and  B 
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I.  Introduction 
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The  velocity  of  the  flow  along  the  plate  is  gradually 
decreasing  due  to  the  wing  installed  on  the  plate.  The 
streamline  is  bent.  The  three-dimensional  turbulent  boundary 
layer  thus  created  will  eventually  separate.  A  separation  region 
will  be  formed  at  the  junction  of  the  wing  and  the  plate. 
Separation  will  lead  to  the  formation  of  a  horseshoe  shaped 
vortex  to  complicate  the  flow  in  this  region.  This  complicated 
phenomenon  often  occurs  on  aircraft  and  other  fluid  equipment, 
affecting  the  performance  of  such  devices.  Flow  separation  and 
its  effect-  on  the  flow  itself  will  affect  the  behavior  of  the 
turbulent  boundary  layer  in  and  out  of  the  separation  region. 
Nevertheless,  we  know  very  little  about  it. 

In  this  area,  Shabaka  and  Bradshaw^1!  did  some  experimental 
work  on  the  decay  of  the  horseshoe  vortex  at  low  speeds. 
McMahon^]  measured  the  velocity  and  turbulence  distribution  of 
the  junction  at  low  speeds.  The  focus  of  this  work  is  placed  on 
the  effect  of  separation  on  the  turbulent  boundary  layer  in  and 
out  of  the  separation  region.  The  barrier  method  used  by  other 
people  are  plates  with  an  elliptical  leading  edge.  Moreover,  the 
point  of  measurement  is  located  far  away  from  the  leading  edge 
and  the  local  pressure  gradient  is  zero.  Hence,  the  behavior  of 
turbulence  measured  does  not  include  any  pressure  gradient 
effect.  In  this  work,  the  straight  wing  used  has  a  thickness 
ratio  of  25%  (See  Figure  1) .  The  purpose  is  to  enhance  the 
strength  of  the  separation  vortex  to  intensify  its  effect  on  the 


turbulent  behavior.  Furthermore,  due  to  the  fact  that  a  straight 
wing  is  used,  the  work  can  more  truthfully  reflect  many  practical 
situations. 

In  this  work,  two  measuring  cross-sections  were  chosen  at 
low  speeds.  One  is  at  18.75%  chord  length.  It  is  135mm  from  the 
leading  edge  of  w:ng  and  is  perpendicular  to  the  chord  (See 
Figure  1) .  It  is  relatively  far  away  from  the  leading  edge  and 
the  upper  and  lower  surfaces  are  parallel.  Thus,  the  pressure 
gradient  is  very  small.  The  other  is  on  the  line  45°  from  the 
chord  near  the  leading  edge.  The  pressure  gradient  and  curvature 
of  streamline  are  relatively  large  over  there.  Five  measuring 
stations  were  selected  on  each  cross-section,  to  be  located  in 
and  out  of  the  separation  region.  At  each  station,  the  average 
velocity  and  the  distribution  of  the  six  turbulent  stress 
components  of  the  turbulent  boundary  layer  were  measured.  At 
subsonic  speeds,  a  station  was  selected  in  and  out  of  the 
separation  region  (See  Figure  2)  to  determine  the  average 
velocity  and  turbulent  stress. 
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Figure  1.  Installation  of  the  Vertical  Wing  Model  and  Positions 
of  Boundary  Layer  Measuring  Stations  at  Low  Speeds. 

1.  vertical  wing 

2.  measuring  station 

3.  plate 

4.  separation  line 

5.  metal  trip  wire 

6.  plate  leading  edge 


Figure  2.  Wing  Section  of  Subsonic  Model  and  Location  of 
Boundary  Layer  Measuring  Stations 

1.  measuring  station 

2.  separation  line 

3.  subsonic  mode  wing  section 


II.  Experimental  Apparatus 
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Wind  Tunnel 

Low  speed  experiments  were  carried  out  in  the  92x92cm 
section  D-2  wind  tunnel  at  Beijing  Institute  of  Aeronautics  and 
Astronautics.  When  the  wind  speed  is  16m/sec,  the  turbulence  in 
the  center  of  the  flow  field  is  0.5%.  Non-uniformity  of  velocity 
pressure  distribution  is  less  than  0.7%.  Deviation  of  flow 
direction  is  less  than  0.5°.  The  Reynolds  number  per  meter  is 
approximately  1x10^. 

Subsonic  experiments  were  done  in  the  G-3  wind  tunnel  at 
Beijing  Institute  of  Aeronautics  and  Astronautics.  This  is  a 
blow  down  wind  tunnel.  Its  subsonic  and  transonic  speed 
experimental  section  is  53. 5x37. 5cm.  The  left  and  right  walls 
are  fixed.  The  upper  and  lower  walls  are  open.  At  M=0.4,  Re 
number  per  meter  during  the  experiment  is  approximately  lxlO7. 

The  turbulence  is«l%.  The  non-uniformity  of  M  distribution  is 


AM 

M 


<0.004 


Vertical  Wing  Model  and  Measuring  Plate 
The  head  of  the  vertical  wing  on  the  measuring  plate  is 
semi-circular.  It  is  immediately  followed  by  a  straight  section 
from  12.5%  to  25%  chord  length.  Afterward,  the  thickness 
gradually  decreases  until  the  trailing  edge.  The  plate  used  to 
make  the  measurements  in  low  speed  experiments  was  horizontally 
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installed  at  22cm  away  from  the  lower  wall  of  the  wind  tunnel. 

The  top  of  the  vertical  wing  model  is  right  next  to  the  upper 
wall  of  the  wind  tunnel.  The  bottom  is  immediately  next  to  the 
measuring  plate.  The  span  of  the  subsonic  model  is  equal  to  the 
width  of  the  experimental  section.  The  model  was  sandwiched 
tightly  between  two  side  walls.  One  of  the  side  walls  was  used 
as  the  plate  to  measure  the  behavior  of  its  boundary  layer.  When 
low  and  subsonic  speed  experiments  were  conducted,  a  0.5mm  metal 
trip  wire  was  placed  near  the  leading  edge  of  the  plate  to  ensure 
that  a  fully  developed  turbulent  boundary  layer  exists  on  the 
plate. 


Measuring  Instruments 

The  average  velocity  and  turbulence  of  the  boundary  layer 
were  measured  by  an  "x"-shaped  dual  hot  wire  probe.  Using  a 
turning  mechanism,  it  was  possible  to  move  the  probe  up  and  down 
perpendicular  to  the  plate  in  the  process.  The  probe  could  also 
be  moved  around  the  point  of  measurement.  This  mechanism  is 
driven  by  a  step  motor  which  can  be  operated  outside  the  wind 
tunnel  by  using  a  controller.  The  minimum  displacement 
controllable  is  0.01mm.  The  minimum  angle  of  rotation  is  0.1°. 

Other  instruments  required  in  the  measurements  include:  an 
anemometer  and  a  barometer  to  measure  wind  speed  and  dynamic 
pressure,  thermocouple  thermometers  to  measure  total  temperature 
of  the  flow  and  wall  temperature,  and  a  Model  TSI  1050  constant 
temperature  hot  wire  anemometer  and  its  accessories,  such  as 
linearizer,  correlator  and  RMS  voltmeter  for  hot  wire 
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measurements.  The  experimental  data  were  read  on  the  RMS 
voltmeter  or  sent  into  a  data  acquisition  system  (Model  XJ-200 
scanner)  which  was  connected  to  a  computer.  In  high  speed 
experiments,  the  data  acquisition  system  can  also  gather  the 
temperature  of  the  stable  section  and  the  wall  temperature  of  the 
experimental  section  with  thermocouples  as  well  as  total  pressure 
of  the  stable  section,  the  static  pressure  on  the  wall  in  the 
experimental  section  and  the  static  pressure  at  the  measuring 
point.  The  data  can  then  be  processed  by  a  computer  (DJS-130) 
before  the  final  results  are  made  available. 

III.  Experimental  Method 

In  order  to  measure  the  three  velocity  components  and  six 

turbulent  stress  components  at  a  point  in  space  with  a  hot  wire 

probe,  it  is  necessary  to  determine  the  orientation  and  position 

of  the  hot  wire.  To  this  end,  three  coordinate  systems  are 

introduced  (See  Figure  3) :  the  wind  tunnel  coordinate  system  x- 

y-z  which  is  fixed  on  the  plate,  the  local  wind  direction 

coordinate  system  s-y-n,  and  the  coordinate  system  T-N-B  which  is 

fixed  on  the  hot  wire.  On  the  plate,  the  x-axis  points  at  the 

incoming  flow  direction;  s-axis  points  at  the  projection  of  the 

local  current  on  the  plate.  The  y-axis  is  perpendicular  to  the 

plate,  pointing  upward.  T-axis  is  along  the  direction  of  the  hot 

wire.  N-axis  is  perpendicular  to  the  hot  wire  in  the  plane 
formed  by  the  T-axis  and  by  the  hot  wire  rotating  axis.  The  included  angle 
between  the  N-axis  and  y-axis  is  a.  The  included  angle  between  s-axis 
and  the  hot  wire  projected  on  the  s-n  plane  is  $ .  The  included  angle 


between  s-axis  and  x-axis  is  0.  Once  the  distribution  of  0  along 
y-direction  is  known,  it  is  possible  to  determine  the  s-axis 
direction  along  y.  In  order  to  measure  0  ,  an  "x"-shaped  hot  wire 
probe  with  a=0  should  be  used.  The  line  equally  dividing  the 
angle  between  the  two  hot  wires  should  coincide  with  the  x-axis. 
The  value  of  ^  can  be  calculated  according  to  the  following: 


where  Ec  A and  Ecg  are  the  linearized  voltage  signals  from  the  two 
hot  wires.  By  choosing  the  proper  4>  and  a,  components  of  average 
velocity  and  turbulence  of  each  point  in  the  T-N-B  coordinate 
system  along  y-direction  in  the  boundary  layer  can  be  obtained. 
Based  on  coordinate  transformation,  their  components  in  the  s-y-n 
coordinate  system  can  then  be  obtained.  Finally,  they  are 
transformed  to  the  x-y-z  coordinate  system.  Thus,  the  average 
velocity  and  turbulence  at  any  point  normal  to  the  plate  in  the 
wind  tunnel  coordinate  system  can  be  obtained. 

Specifically,  experimental  data  are  processed  as  follows. 

The  effective  velocity  to  cool  the  hot  wire  is  introduced  as 
Oeff.  Then,  the  linearized  voltage  output  is  linear  in  relation 
with  varying  p Deff?  i.e.. 


E,= 


pU  .f  i 

K 


(2) 


where  K  is  a  proportionality  constant  which  can  be  determined  by 
calibrating  the  hot  wire.  Obviously,  the  cooling  speed  Ueff 


ought  to  be  a  function  of  UN,  UT  and  UB.  Here,  let  us  introduce 
the  expression  used  by  Jorgensen ^1. 

£/?+*«  £/’)■/>  (3) 

where  UT,  CN  and  UB  are  velocity  components  along  T-N-B  axes, 
respectively.  Coefficients  c  and  h  still  must  be  determined  by 
calibrating  the  hot  wire.  If  the  orientation  of  the  hot  wire, a 
and  is  known,  it  is  very  easy  to  derive  the  relation  between 
the  velocity  components  UT,  UN,  UB  and  their  corresponding 
velocity  components  us+u's,  Uy+Uy,  u^  on  the  s-y-n  coordinate 
system  (where  us  and  Uy  are  average  values,  and  u's  and  u^  are 
fluctuations) .  Based  on  equation  (2) ,  Ec  can  be  expressed  as 

Ec=Ec[p(us+Ug),  P(uy+Uy),  P\i'ni  *,  a]  (4) 

and 

E.^E'+e,  (5) 

where  Ec  and  ec  are  the  mean  and  fluctuating  voltage  output  from 
the  linearizer,  respectively.  Let  us  take  the  time  average  of 
both  sides  of  equation  (4)  and  also  take  the  time  average  of  the 
square  of  both  sides  of  equation  (4) ,  we  get  (after  omitting 
higher  order  terms)  the  following: 


where  the  average  flow  density  p"  can  be  obtained  based  on  ~p=P/RT 
by  measuring  the  static  pressure  P  and  mean  temperature  T. 
Parameters  A,  B,  C,  D,  E  and  F  are  functions  of  ^  and  a.  All 
velocity  fluctuating  terms  are  mass  weighted  to  the  following 
definitional  : 


u.l  =  (pu',)*/ p*  ,(u\  u,),  =  p*  u ,  u,  / p' 

Similar  definitions  apply  to  terms  with  different  subscripts 
(s,y,n)  . 


Figure  3.  Wind  Tunnel,  wind  Direction  and  Hot  Wire  Coordinate 
Systems 
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1.  projection  of  hot  wire  on  s-n  plane 

2.  hot  wire 

3.  axis  of  rotation 


There  are  eight  unknowns  in  equation  (6)  and  (7) .  Two  of 
them  are  average  velocity  components:  us  and  uy  (un=0) .  Six  of 
them  are  velocity  fluctuation  related  terms  which  take  density 
variation  into  account  (to  reflect  six  turbulent  stress  terms) : 

\ »  fQ  ( u'.u ,  TmTOT.  iu'.u',). 

If  different  hot  wires  are  used  and  proper  and  ex  are  chosen 
based  on  the  feasibility  to  be  installed  in  the  wind  tunnel, 
equation  (6)  and  (7)  can  be  simplified  by  various  degrees.  Thus, 
experiments  can  be  performed  with  respect  to  certain  and  a  to 
solve  equation  (6)  and  (7)  by  using  the  values  of  Ec  and  e£, 
measured.  The  eight  unknowns  can  thus  be  determined.  Through 
coordinate  transformation,  the  average  velocity  and  turbulence 
along  the  wind  tunnel  coordinate  system  can  be  obtained: 


|n  u,\  ,  u,'  ,  u,\  .  (u'.u',). 


It  can  be  proven  that,  under  similar  conditions  to  derive 
equation  (6)  and  (7)  ,  it  is  possible  to  get  u'yp  and  u^  from 

u'x2,  uy2  and  u'z2.  The  correlation  is 


and 


where  cr=(y-l)M2  and  the  specific  heat  V= 1.4.  As  for  the 
correlating  experimental  data,  with  the  exception  of  shear  at 
M=»0 . 4 ,  the  rest  are  expressed  by  conventional  turbulent  stress 


terms. 


IV.  Analysis  of  Experimental  Results 


Based  on  the  experimental  results,  near  the  separation 
region,  the  average  speed  and  turbulent  stress  are  both  affected 
by  the  main  vortex  of  separation.  It  is  more  apparent  near  the 
main  vortex.  In  order  to  clearly  illustrate  this  effect,  the 
experimental  curve  was  organized  by  x,y,z  to  directly  reflect  the 
actual  position  of  the  separation  vortex. 

(1)  First,  let  us  discuss  the  results  measured  on  the 
x=135mm  cross-section  away  from  the  leading  edge  of  the  wing. 

The  pressure  gradient  of  this  cross-section  is  dp/dx«0. 

Moreover,  the  streamline  is  very  straight.  Therefore  the 
experimental  data  obtained  is  not  influenced  by  pressure  gradient 
and  streamline  curvature.  In  this  case,  some  of  the 
characteristics  of  a  uniform  flow  can  be  visualized  from  the 
velocity  vector  distribution  projected  on  the  cross-section  (See 
Figure  4a) .  Based  on  the  core  position  of  the  main  separation 
vortex  (vertical:  y=13mm,  spanwise:  z=50mm)  and  the  induced 
velocity  caused  by  the  separation  vortex  shown  in  this  figure,  we 
can  see  that  the  induced  velocity  increment  is  very  large  in  the 
z-direction  near  the  plate.  This  induced  velocity  causes  high 
speed  flow  away  from  the  wing  to  be  rolled  near  the  vertical 
wing.  The  boundary  layer  thickness  is  thus  reduced.  Also,  low 
speed  current  near  the  wing  is  rolled  away,  causing  the  boundary 

layer  to  become,  thicker  there,  and,  furthermore,  the  distribution 

in  the  y  direction  of  ux  causes  fluctuation  in  u*  distribution  near  the 
main  vortex  (See  Figure  4b). 
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(a)  (b) 

Figure  4 

(a)  velocity  component  distribution  on  x=135mm  cross-section 

(b)  ux/u0  y  curve  at  various  stations  at  x=135mm 

1.  vortex  center 

2.  plate 

3.  plate 


Figure  5.  A u' ^ /  u^y  and  -uxu'/u^~y  Curves  at  Various 

•  Stations  at  x=135mm  y 

1.  plate 

2.  plate 

3.  vortex  center 

4.  plate 

5.  plate 

6.  vortex  center 


As  for  quantities  expressing  the  variation  of  turbulent 
stress,  n/u^/u^vs.  y  and  -u'xUy/u£  vs.  y  are  plotted  as  examples 
in  Figure  5  (See  references  [5]  and  [6]  for  other  turbulent 
related  curves).  In  the  figure,  results  of  the  plate  without  the 
wing  are  also  presented  for  comparison  in  order  to  see  the  effect 
of  separation  vortex  on  turbulent  stress  distribution.  Near  the 
main  vortex  (such  as  shown  by  the  two  curves  z=40mm  and  60mm), 
the  curve  shows  a  peak  near  the  core  (y=13mm) .  Near  the  vertical 
wing  (such  as  shown  by  z=10mm) ,  turbulence  stress  decreases 
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significantly.  As  compared  to  the  results  in  reference  [2] ,  the  /45 
peak  turbulent  stress  obtained  in  this  work  is  much  higher. 

However,  the  turbulent  stress  near  the  vertical  wing  is  much 
smaller  compared  to  that  in  reference  [2] .  This  should  be 
related  to  the  fact  that  the  thickness  of  the  model  used  in  this 
experiment  is  larger  so  that  the  separation  vortex  strength  is 
consequently  higher.  In  conclusion,  the  induction  of  separation 
vortex  causes  the  turbulence  to  increase  near  the  main  vortex. 
Consequently,  far  away  from  the  main  vortex,  especially  near  the 
wing,  the  turbulence  decreases.  This  phenomenon  also  appears  in 


vV,’/u.,y^K  and 


vs.  y  curves  (See  reference  [5]). 

(2)  Next,  let  us  discuss  the  results  obtained  near  the 
leading  edge  in  a  cross-section  45°  from  the  chord.  In  this 
cross-section,  the  pressure  gradient  and  streamline  curvature  are 
high.  Hence,  the  results  are  significantly  affected  by  pressure 
gradient  and  streamline  curvature. 

With  regard  to  the  average  boundary  layer  flow  velocity 
there,  Ujj/u,,,  vs.  y  and  uz/Uoo  vs.  y  curves  are  selected  for 
interpretation.  Figure  6  shows  the  distribution  of  the  average 
velocity  component  along  y-direction  in  stations  at  various 
distances  5,.  In  this  cross-section,  it  is  obvious  that  the 
average  velocity  near  the  wing  should  be  higher  than  that  away 
from  the  wing.  This  result  is  apparently  reflected  by  the  group 
of  curves  of  ux/uM  and  uz/uoo  • — •  y.  When  y  is  small,  (near  the 
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plate)  ,  Uz/Uoo  curves  vary  in  a  wave  form  which  is  a  result  of  the 
induction  effect  of  the  separation  vortex.  Based  on  the  peak  of 
the  curve/  the  vortex  core  position  is  approximately  located  at  y« 
4mm.  Similarly/  the  induction  of  separation  vortex  also  causes 
ux/u*  y  curves  to  fluctuate  near  the  vortex  (such  as  shown  by  the 
curve  at  l=40mm) .  The  position  of  the  vortex  center  should  be 
located  between  1=40  and  60mm.  The  distribution  of  experimental 
points  on  the  u^j/u,,,,  ~  y  curve  is  similar  to  that  obtained  with 
stations  on  the  x=135  cross-section.  Due  to  the  effect  of  large 
pressure  gradient  and  streamline  curvature, it  is  even  irore  different 
compared  to  the  plate  curve.  Away  from  the  wing  (such  as 
shown  by  curves  at  1=120  and  160mm) ,  the  disagreement  between 
those  curves  and  the  plate  curve  indicates  that  the  velocity 
model  is  affected  by  the  negative  pressure  gradient. 


c 

Figure  6.  ux/i^~y  and  uzi^~y  Curves  Near  the  Leading  Edge 
in  a  Cross-section  45°  From  the  Chord. 

1.  plate 

2.  plate 

3 .  symbol 

4.  vortex  center 

5.  vortex  center 


As  for  quantities  related  to  turbulent  stress,  distribution 
curves  of  \J uxVu»  and  -u'xUy/u|  were  selected  and  plotted  in 
Figure  7.  In  the  same  figure,  the  experimental  curve  for  the 
plate  is  also  shown  to  compare  the  effect  of  the  wing  on 
turbulent  flow  behavior.  From  the  figure  we  can  see  that  the 
curves  show  peaks  near  the  vortex  center.  In  addition,  the  peaks 
are  much  larger  than  those  obtained  at  the  x=135mm  cross-section. 
Based  on  this  fact,  we  can  conclude  that  the  vortex  has  a 
stronger  effect  on  shear  stress.  The  difference  between  curves 
near  the  wing  (such  as  l=15mm)  and  the  experiment  plate  curve  is 
not  as  large  as  that  at  x^lSSmm.  These  phenomena  show  that  the 
area  affected  by  vortex  is  smaller  near  the  leading  edge.  The 
range  affected  continues  to  expand,  however,  as  we  go  further 
downstream.  Hence,  negative  pressure  gradient  of  the  circulating 
flow  has  a  great  impact. 


Figure  7. 
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v ux2/Uco~y  and  -u'  u'  /u£  y  Obtained  at 
45°  from  the  Chord  Wear  the  Leading  Edge 


1.  vortex  center 

2.  plate 

3.  plate 

4.  vortex  center 

5.  plate 

6.  plate 

(3)  In  order  to  understand  the  changes  in  turbulent  stress 
behavior  in  and  out  of  the  separation  region  at  subsonic  speeds, 
a  station  outside  the  separation  region  at  x=32,  y=15mm,  and 
another  station  in  the  separation  region  at  x=17,  y=15mm  (See 
Figure  2)  were  chosen  to  make  measurements  at  subsonic  speeds. 

At  M=0.4,  \/ ux2/ue>  \J uz2/ue  and  -(u'xu'z)p/ue2  vs.  y/§  curves  were 
plotted  as  shown  in  Figure  8.  Based  on  this  figure,  the 
separation-  vortex  is  approximately  located  at  y/S«0.1.  Due  to 
the  effect  of  the  separation  vortex,  turbulent  stress  increases 
rapidly  and  peaks  there.  This  phenomena  is  basically  the  same  as 
that  at  low  speeds.  Similarily,  its  effect  on  shear  stress  is 
higher  quantitatively  as  compared  to  that  on  positive  stress. 


Figure  8.  Distribution  of  Turbulent  Stress  Related  Parameters 
at  Two  Stations  In  and  Out  of  Separation  Region  at 
M  -0.4. 


V.  Conclusions 


An  "x"-shaped  twin  hot  wire  probe  can  be  used  to  measure  the 
boundary  layer  of  separation  flow  at  low  and  subsonic  speeds 
through  the  proper  choice  of  orientation  angles  ^  and  a.  Related 
characteristic  data  can  be  obtained  by  an  appropriate  data 
processing  technique.  The  experimental  results  show  that  the 
horseshoe  vortex  caused  by  flow  separation  is  an  important  factor 
affecting  the  three-dimensional  turbulent  boundary  layer  behavior 
in  and  out  of  the  separation  region. 

(1)  Effect  of  Separation  Vortex  on  Mean  Velocity  Distribution 

The  induction  effect  of  the  separation  vortex  causes  the 

boundary  layer  thickness  to  decrease. inside  the  separation 
region,  and  it  also  causes  the  boundary  layer  thickness  to 
increase  out  of  the  separation  region.  Near  the  separation 
vortex,  the  average  velocity  shows  fluctuations.  uz  increases 
rapidly  near  the  plate.  The  effect  of  pressure  gradient  on 
average  velocity  is  similar  to  that  on  the  velocity  distribution 
of  the  boundary  layer  of  the  plate. 

(2)  Effect  of  Separation  Vortex  on  Turbulent  Stress 
Distribution. 

The  separation  vortex  decreases  the  turbulent  stress  near  the 
wing.  It  increases  the  turbulent  stress  near  the  vortex  center 
and  a  peak  exists.  The  effect  of  separation  vortex  on  shear 
stress  is  higher  than  that  on  positive  stress.  Near  the  leading 
edge  of  the  wing,  the  separation  vortex  effect  is  relatively 
concentrated  around  the  vortex  center.  As  we  go  further 
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downstream,  the  area  affected  by  the  vortex  expands  larger.  The 
strength  decreases  and  the  peak  drops.  In  this  case,  the 
negative  pressure  gradient  has  some  effect. 

(jj  At  subsonic  speeds,  the  effect  of  separation  vortex  on 
turbulent  stress  distribution  is  the  same  as  that  at  low  speeds. 
The  effect  of  the  vortex  on  shear  stress  is  also  higher  than  that 
on  positive  stress. 

Comrade  Ma  Shulin  of  Beijing  Institute  of  Aeronautics  and 
Astronautics  also  worked  on  this  project.  Senior  Students  Wang 
Yun  and  Lu  Gan  of  Aerodynamics  Major  at  Beijing  Institute  of 
Aeronautics  and  Astronautics  of  the  Class  of  1983  worked  on  the 
low  speed  experiments.  Graduate  student,  Wang  Chaoan  worked  on 
the  high  speed  experiments. 
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Effect  of  Winglet  on  Spatial  Vortex  of  a  Solid  in  Rotation  at  /49 

High  Angle  of  Attack 
Wang  Zixing  and  Wu  Genxing 
(Nanjing  Aeronautical  Institute) 

Abstract 

This  paper  presents  some  experimental  results  on  the 
mechanism  of  the  effect  of  the  winglet  on  the  spatial  vortex  of  a 
rotating  body.  The  results  show  that,  with  the  winglet,  for  a 
rotating  body  without  skidding,  the  position  of  the  spatial 
vortex  is  lowered  because  its  strength  is  reduced.  Hence,  the 
asymmetric  problem  at  high  angle  of  attack  is  resolved  somewhat 
by  various  degrees. 


I.  Introduction 

Lateral  force  which  occurs  at  high  angle  of  attack  without 
side  skidding  is  primarily  due  to  asymmetric  vortices.  In  order 
to  eliminate  and  delay  the  appearance  of  this  undesirable  lateral 
force,  many  researchers  performed  force  measurements  in  wind 
tunnels  and  presented  several  solutions.  For  example,  a  strake 
in  the  front  was  used  ^-6] #  However,  there  is  relatively  little 
understanding  on  the  vortex  system.  On  the  basis  of  reference 
[2],  a  fluorescent  mini-tuft  technique^  was  used  in  this  work 
to  clarify  the  spatial  vortex  and  winglet  vortex  of  a  rotating 
solid  in  order  to  understand  the  role  of  the  winglet  in 
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eliminating  asymmetric  vortices. 

Five  different  winglet  types  and  arrangements  were  used  to 
study  the  spatial  vortex  of  a  sharp  rotating  body  at  high  angle 
of  attack  in  the  absence  of  side  skid  in  a  low  speed  wind  tunnel. 
The  results  show  that  after  being  equipped  with  various  winglets,  the 
asymmetric  vortices  of  the  body  of  rotation  at  high  attack  angles 
can  be  suppressed  to  some  extent.  The  best  result  was  obtained 
with  a  small  strake  at  high  sweepback  angles. 

We  also  used  the  conventional  oil  flow  pattern  method  to 
observe  the  changes  of  separation  lines  after  winglets  were 
attached  to  some  models. 

II.  Experimental  Apparatus,  Model  and  Method 

1.  Experiments  were  carried  out  in  a  0.75m  opening  wind 
tunnel.  The  wind  speed  was  30  m/sec.  When  the  model  diameter 
was  used  as  the  reference  length,  the  Reynolds  number  Re=1.2xl05. 

2.  The  position  of  the  spatial  vortex  was  determined  with 
the  assistance  of  a  three  coordinate  displacement  mechanism  which 
was  equipped  with  fluorescent  tuft  probes.  Details  of 
fluorescent  tuft  visualization  technique  and  spatial  vortex 
measurement  are  shown  in  reference  [1],  The  coordinate  of  the 
vortex  trajectory  was  rendered  dimensionless  by  using  the  length 
of  the  body  of  rotation  as  a  reference. 

3.  The  fuselage  size  of  all  five  different  models  remains 
the  same.  The  total  aspect  ratio  is  7.  The  aspect  ratio  at  the 
tip  is  3.5.  The  rear  part  of  the  model  is  cylindrical.  The 
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front  of  the  model  is  pointed,  made  of  a  circular  arc.  The 


dimension  and  position  of  the  winglet  is  shown  in  the  figure 


below.  The  wing  is  Model  NACA  64A  004.  The  installation  angle 


of  the  wing  surface  (with  the  exception  of  Model  5)  is  0°. 


4.  Experimentation  began  from  15°  and  ended  at  60°. 


Manuscript  received  on  May  28,  1984.  Revised  manuscript  received 
on  September  5. 
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1.  Model  1 


2.  Model  2 


3.  Model  3 


4.  Model  4.  Winglet  installed  at  0( 


5.  Model  5.  Winglet  installed  at  50°.  The  middle  of 


the  wing  root  chord  is  used  as  the  axis  of  rotation, 
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Figure  1 

1.  body  vortex  without  winglet 

2.  body  vortex  with  winglet 

3.  tail  vortex  of  winglet 

4.  gradually  diffusing 

5.  asymmetric  body  vortices 

6.  body  vortex  without  winglet 

7.  tail  vortex  of  winglet 

8.  body  vortex  with  winglet 

9.  (Note:  Body  vortices  shown  without  winglet  are  side 

views. ) 

10.  Model  l,a*  22.5° 

11.  Model  1,  a=  45° 

12.  Body  vortex  without  winglet 

13.  body  vortex  with  winglet 

14.  tail  vortex  of  winglet 

15.  asymmetric  vortices  without  winglet 

16.  Model  1,  a  -  50° 
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III.  Results  and  Analysis 

1.  Small  Strake  Installed  at  the  Front  of  a  Slender  Body 
(Model  1) 

Figure  1  shows  how  the  spatial  vortices  of  a  slender  body 
vary  with  angle  of  attack  with  and  without  a  small  strake.  With  / 
a  small  strake,  when  the  attack  angle  is  not  very  large  (such  as 
a  =10°) .  because  the  tail  vortex  of  the  winglet  is  very  close  to 
the  axis  of  the  slender  body,  its  induction  effect  delays  the 
separation  of  the  slender  body.  Consequently,  spatial  vortex  did 
not  form.  •  When  a =22. 5°,  due  to  the  presence  of  winglet,  body 
vortex  could  only  be  formed  behind  the  wing.  Its  strength  is 
less  than  that  of  a  wingless  slender  body.  The  position  is  also 
lower. 

As  the  angle  of  attack  increases  (a>30°),  the  wing  vortex 
diffuses  gradually.  The  position  and  strength  of  body  vortices 
are  obviously  suppressed  by  the  winglet  so  that  they  cannot 
develop  into  an  asymmetric  vortex  system.  The  tail  vortex  of  the 
winglet  did  not  roll  into  the  body  vortex. 


Figure  2.  Separation  Lines  With  and  Without  Winglet 

1.  separation  line  without  winglet 

2.  separation  line  with  winglet 

Figure  2  was  drawn  based  on  oil  flow  pattern  pictures.  From 
the  figure  we  can  see  that  due  to  the  presence  of  the  winglet  the 
position  of  separation  is  obviously  lowered.  The  development  of 
the  body  vortex  was  delayed. 

On  the  rear  part  of  the  model,  in  addition  to  a  pair  of 
primary  vortices,  there  are  a  pair  of  secondary  vortices  in 
opposite  directions  and  a  pair  of  relatively  weak  secondary 
vortices.  In  order  to  avoid  confusion  due  to  too  many  lines,  the 
opposite  secondary  and  secondary  vortices  are  omitted.  Only  the 
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Figure  3 

(a)  Model  2,  «=  35°  (b)  Model  2,  <x-  45°  (c)  Model  2,  50° 

•* 

1.  asymmetric  pair  of  vortices  without  winglet 

2.  gradual  breakdown 

2.  Small  Strake  Installed  in  Mid-Section  of  Slender  Body 
(Model  2). 

Based  on  Figure  3  we  can  see  that  when  the  attack  angle  is 
high  (such  as  <s>22.5°),  there  are  two  strong  front  body  vortices 
produced  by  the  head  of  the  body  of  rotation.  In  addition,  due  to 
the  presence  of  the  winglet,  a  very  strong  wing  vortex  exists  in 
the  mid-section.  The  body  vortices  are  lowered  by  the  wing 
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vortex.  Because  of  the  presence  of  the  wing,  the  vortex 
generated  by  the  slender  body  at  the  wing  is  not  rolled  into  the 
front  body  vortex.  Instead,  it  is  rolled  into  a  new  body  vortex 
behind  the  wing.  Due  to  the  influence  of  the  wing  vortex,  the 
position  and  strength  of  the  body  vortex  are  suppressed 
significantly.  Thus,  the  asymmetry  of  body  vortex  is  alleviated. 
When  the  attack  angle  is  greater  than  50°,  both  wing  and  body 
vortices  are  breaking  down.  Hence,  there  are  no  asymmetric 
vortices.  When  a-  60°,  wing  vortex  totally  breaks  down  and  body 
vortices  exist  only  at  the  head.  Due  to  the  effect  of  the  low 
pressure  zone  caused  by  the  breakdown  of  the  tail  vortex  of  the 
winglet,  Lt  does  not  swing  upward  as  in  the  case  where  there  is 
no  winglet. 

Figure  4  shows  the  effect  of  the  winglet  on  the  separation 
line  and  on  body  vortices. 


i.  '■* 
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Figure  4 

(a)  Model  2,  a  =  45°  (b)  Model  2,  a =  45°  Top  View 


1.  separation  line 


3.  Trapezoidal  Strake  With  Small  Span  Installed  in 
Mid-Section  of  Slender  Body 

The  situation  with  Model  (3)  is  similar  to  that  with  Model 
(2).  The  winglet  causes  the  body  vortex  to  drop.  When  <*=  30°, 
body  vortices  are  merged.  As  a  increases  further,  the  wing 
vortex  breaks  down.  Body  vortices  also  break  because  of  entering 
the  breakdown  zone.  As  compared  to  a  triangular  winglet,  the 
separation  breakdown  of  this  type  of  winglet  occurs  at  a  smaller 
attack  angle.  Hence,  it  is  less  effective  in  eliminating  the 
asymmetric  vortices.  However,  it  will  not  cause  any  asymmetry  at 
all. 

4.  Small  Trapezoidal  Winglet  Installed  in  the  Middle  of 
the  Slender  Body  (Model  4) 

The  wing  vortex  produced  by  the  winglet  also  significantly 
affects  the  body  vortex  to  deflect  it  lower  and  to  weaken  it 
(because  it  cannot  superimpose  onto  itself  in  an  area  where  the 
winglet  is).  When  a=22.5°,  the  body  vortex  is  induced  by  the 
wing  vortex  and  disappears.  The  wing  vortex  moves  downstream. 
When  a«30°,  the  wing  vortex  breaks  down.  The  right  and  left  body 
vortices  are  symmetric.  They  also  gradually  begin  to  break  down 
(See  Figure  5) . 


Figure  5 


Model  4,  a =30°,  Wing  vortex  already  broken  down, 
symmetric  body  vortices 


5.  Small  Trapezoidal  Wing  Installed  in  the  Middle  of  the 
Slender  Body  at  an  Angle  (Model  5) 

Because  it  was  found  experimentally  that  separation  also 
began  at  a =30°  with  Model  4,  in  this  test  the  winglet  was 
installed  at  -30°.  When  a=35°/-~40°,  due  to  the  fact  that  the 
attack  angle  of  the  winglet  is  only  5°/^10°,  the  wing  vortex  was 
relatively  weak.  It  had  little  effect  on  the  body  vortex.  At  a 
*45°,  the  wing  vortex  became  stronger.  The  body  vortex  was 
induced  downward.  As  a  increased  further,  the  wing  vortex 
intensified  more.  The  body  vortex  disappeared.  Based  on  this 
result,  as  long  as  the  actual  attack  angle  of  the  winglet  is  not 
greater  than  30°,  it  can  effectively  suppress  the  asymmetric 
separation  vortex  sy3tem  of  a  slender  body  at  large  attack  angles 
in  the  absence  of  side  skidding. 


IV.  Conclusions 
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This  paper  introduced  a  study  of  separation  vortex  systems 
of  five  types  of  slender  bodies  with  winglets  in  a  low  speed  wind 
tunnel  using  a  fluorescent  tuft  technique.  A  intuitive 
understanding  on  the  effective  suppression  of  the  development  of 
asymmetric  vortex  systems  by  the  winglet  was  obtained. 

1.  A  winglet,  in  addition  to  being  capable  of  suppressing 
the  development  of  body  vortex  (partly  because  there  is  no  input 
and  partly  because  the  separation  position  is  lowered)  near 
itself  (shaded  area  towards  downstream) ,  when  its  spatial  vortex 
has  considerable  strength,  can  effectively  attract  the  separation 
vortex  of  the  slender  body  so  that  it  will  not  be  deflected 
upward  to  develop  into  an  asymmetric  vortex  system.  Thus,  the 
undesirable  lateral  force  can  be  avoided  in  the  case  side 
skidding  is  absent. 

2.  Separation  occurs  when  the  attack  angle  of  the  winglet 
is  too  high.  Although  this  is  unfavorable  for  drag,  the  spatial 
vortex  is  sucked  into  the  separation  zone.  Thus,  it  is  not 
possible  to  develop  into  an  asymmetric  vortex  system. 

3.  Any  bulge  on  either  side  of  the  slender  body  can  reduce, 
delay  or  eliminate  the  formation  of  an  asymmetric  vortex  system 
because  its  tail  vortex  can,  to  some  extent,  suppress  and  attract 
the  separation  vortex  of  the  slender  body. 
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On  the  Investigation  of  Separated  Flow  in  the  Base 
Bian  Yinqui  (Institute  of  Mechanics)  and 
Dong  Changquan  (Beijing  Institute  of  Information  and  Control) 

I.  Introduction 

Separated  flow  in  the  base  is  an  important  topic  in 
aerodynamics.  In  engineering  applications,  the  primary  work  is 
to  calculate  the  pressure  and  heat  flow  in  the  base  in  order  to 
provide  a  qualitative  as  well  as  quantitative  description  of  the 
flow  field  structure.  Due  to  the  complexity  of  the  base  flow 
field  itself,  there  are  not  enough  experimental  data  especially 
in  hypersonic  range.  Methods  for  theoretical  analysis  have  not 
yet  been  perfected.  In  particular,  the  study  of  heat  transfer  in 
the  base  is  far  behind  as  compared  to  the  study  of  pressure  in 
the  base.  This  is  because  the  heating  process  cannot  be 
reasonably  treated  easily  under  the  condition  that  the  base  flow 
field  structure  is  not  known.  Some  empirical  methods  are  subject 
to  limitations.  Direct  solution  of  the  complex  N-S  equations  by 
numerical  methods  has  not  yet  reached  the  level  of  providing 
reliable  design  data.  Hence,  designers  are  caught  in  the  middle 
when  choosing  the  base  heat  flow  data.  Based  on  the  nature  of 
the  flow  field,  the  flow  in  the  base  has  many  things  in  common 
with  the  separated  flow.  Similar  analytical  treatments  and 
numerical  methods  can  also  be  referenced  either  way.  Regardless 
whether  it  is  under  the  influence  of  an  external  flow  or  due  to 
geometric  changes  on  the  surface,  the  flow  is  separated  when  a 


high  pressure  area  appears.  In  addition,  a  re-circulation  zone 
is  formed  near  the  wall.  Later,  the  flow  rejoins  the  surface. 
Similarly,  separated  flow  in  the  base  also  includes  an  external 
inviscid  flow,  a  free  shear  layer  and  a  re-circulation  zone.  T! 
free  shear  layer  is  similar  to  the  surface  boundary  layer.  The 
difference  is  that  both  boundaries  of  the  shear  layer  are 
unknown.  Figure  1  shows  a  schematic  diagram  of  the  separated 
flow  in  the  base. 


Figure  1.  Schematic  Diagram  of  Base  Separated  Flow 

1.  boundary  layer 

2.  inviscid  flow  region 

3.  shock  wave 

4.  expansion  wave 

5.  free  shear  layer 

6.  boundary  streamline 

7.  trailing  shock  wave 

8.  turning  region 

9.  expansion  core 

10.  rear  stationary  point 

11.  re-circulation  region 

12.  base  region 

13.  dull  body 


The  study  of  engine  jet  is  another  base  flow  problem.  The 
special  feature  of  a  jet  flow  field  is  that  it  is  usually  a 
turbulent  flow.  Furthermore,  it  is  a  mixture  of  chemically 
reacting  gases.  A  schematic  diagram  of  the  jet  is  shown  in 
Figure  2.  The  contents  of  the  study  include  chemical  kinetics, 
turbulent  modeling,  and  numerical  simulation. 

A  simplified  analysis  usually  has  the  following  assumptions: 
the  exit  flow  field  is  symmetric  with  respect  to  the  parallel 
axis.  If  the  base  diameter  is  much  larger  than  the  nozzle 
diameter,  the  base  re-circulation  zone  will  not  significantly 
affect  the  entire  jet  structure.  In  addition,  it  is  also  assumed 
that  there  is  no  radial  pressure  gradient.  Furthermore,  some 
complicated  factors  are  also  not  considered;  such  as  the  effect 
of  shock  wave,  thermal  relaxation  of  condensed  particles  and 
velocity  relaxation.  Thus,  the  jet  is  just  as  an  ordinary 
stream.  Equations  with  boundary  layer  conditions  are  used; 
including  mean  shear  layer  equations  using  velocity,  temperature 
and  composition  as  variables.  In  China,  a  few  organizations  have 
started  some  experimental  work.  There  is  very  little  theoretical 
analysis  done.  Work  done  abroad  only  began  in  recent  yearst*”^ . 

Manuscript  received  on  July  16,  1984. 


Generalized  separated  flow  in  the  base  should  also  include  /55 
the  study  of  the  complicated  wake  formed  due  to  the  detachment  of 
the  vortex.  The  flow  around  a  dull  body  has  a  wide  range  of 
applications.  There  are  topics  in  this  area  in  aeronautics, 
aerospace,  ocean  engineering  and  wind  load  of  tall  buildings. 

There  is  no  profound  understanding  on  the  mechanism  of  vortex 
detachment  behind  a  dull  body.  Nevertheless,  experimental  work 
is  way  ahead  of  theoretical  analysis.  Back  in  the  sixteenth 
century,,  before  fluid  dynamics  was  established  as  a  discipline, the 
famous  Italian  painter  and  scientist,  Leonardo  da  Vinci  (1452- 
1519)  drew  pictures  showing  the  re-circulation  and  vortex  of 
flows  toward  a  plate,  around  a  column  and  around  a  triangular 
beam^J't5!.  It  can  be  considered  as  the  earlier  record  noting 
the  observation  and  study  of  vortex  detachment.  Earlier  in  this 
century,  Planck  took  some  photographs  of  the  wake  vortex  motion 
around  a  circular  beam  .  These  are  qualitative  descriptions 
based  on  flow  field  visualization.  There  was  very  little 
theoretical  analysis  in  earlier  days.  As  testing  techniques  were 
developed  and  computers  became  available,  together  with  a  strong 
background  in  applications,  a  great  deal  of  progress  was  made  in 
areas  involving  the  vortex  formation  mechanism  behind  the  dull 
body  and  the  interaction  of  vortices  Let  us  use  the 

flow  around  a  circular  column  as  an  example,  it  is  the  most 
widely  studied  topic  in  fluid  dynamics  experimentally, 
theoretically  and  in  computation.  The  geometry  of  a  two- 
dimensional  circular  column  is  very  simple.  However,  even  in  a 
steady  flow,  the  near  wake  flow  field  varies  in  a  complex  manner 
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in  various  Reynolds  number  ranges.  The  mechanisms  for  vortex 
detachment,  vortex  interaction  and  vortex  streak  formation  are 
not  yet  known.  Figure  3  shows  the  experimental  results  obtained 
over  the  years^^]  .  It  represents  the  stationary  point  pressure 
coefficient  Cpb  vs.  Re. 


Figure  2.  Schematic  Diagram  of  Engine  Jet 

1.  core  area 

2.  mixing  area 

A-B  (Re=l-50) :  The  entire  flow  field  is  a  steady  laminar 

flow. 

B-C:  If  Re >50  and  the  flow  remains  to  be  steady  and 

laminar,  based  on  theoretical  analysis,  Cpb  varies  along  the 
curve  B-C. 

B-D  (Re=50-200) ;  the  actual  situation  does  not  follow  B-C. 
Although  the  flow  field  still  remains  laminar,  the  wake  begins  to 
oscillate  periodically. 

D-E  (Re=200-1500) :  Turbulence  begins  to  appear  and  periodic 
oscillation  moves  downstream. 

E-F  (Re=1500-2xl0^) ;  The  turning  point  moves  forward  along 
the  free  shear  layer.  The  turning  process  develops  into 
turbulent  flow  from  periodic  oscillation. 


Figure  3.  Pressure  Coefficient  at  Stationary  Point  Behind 
Circular  Column 

1.  steady  laminar  flow 

F-G  (Re=5xl05+) :  In  this  region,  the  drag  index  falls 
drastically  due  to  rising  base  pressure.  The  turning  point  is 
very  close  to  the  separation  point.  The  flow  always  starts  with 
laminar  separation.  After  forming  a  bubble  zone  on  the  surface 
due  to  attachment,  it  becomes  turbulent  separation. 

The  pressure  coefficient  plot  is  somewhat  typical.  However, 
it  is  significantly  affected  by  the  turbulence  of  the  main  flow, 
the  roughness  of  the  surface  and  other  factors.  It  is  not 
possible  to  discuss  this  problem  in  depth  here.  Similar  topics 
were  also  discussed  in  the  "Meeting  on  Separation  Flow  and  Vortex 
Motion"  and  reference  [9]  was  introduced.  In  this  paper,  we  will 
briefly  introduce  the  experimental  research  and  theoretical 
analysis  on  the  base  flow  of  flight  vehicles  in  Section  II.  In 
Section  III,  we  will  concentrate  on  solving  N-S  equations  with 


numerical  methods 


II.  Brief  Introduction  to  Experimental  Research  and  Analytical 


Methods  on  Base  Flow 

Most  base  pressure  experiments  are  done  at  low  Mach  numbers 
(M<5) .  There  are  few  high  Mach  number  experiments  and  even  fewer 
of  them  are  done  in  turbulent  flow.  Because  it  is  difficult  for 
ground  simulation  equipment  to  simultaneously  satisfy  the  high 
Mach  number  and  high  Reynolds  number  requirement,  a  great  deal  of 
progress  has  been  made  continuously  in  the  are  of  testing 
techniques.  The  experimental  range  can  be  extended  to  around 
M=30  and  Re=107~108. 

For  slender  cones,  there  are  several  parameters  affecting  the 
base  pressure;  including  laminar  or  turbulent  flow,  Mach  number 
and  Re  number  of  the  incoming  flow,  geometric  parameters  of  the 
cone  (semi-conical  angle  9C  and  dullness  Rfl/Rb  where  RN  is  the 
radius  of  the  head  and  R^  is  the  radius  of  the  base) ,  specific 
wall  temperature  and  the  magnitude  of  the  attack  angle. 

(1)  Effect  of  Re:  The  state  of  the  flow  (laminar  or 
turbulent)  is  an  important  factor  determining  the  base  pressure. 
In  the  case  of  laminar  flow,  the  base  pressure  drops  with 
increasing  Re.  In  the  case  of  turbulent  flow,  it  is  independent 
of  Re.  Nevertheless,  some  individual  experiments  found  that  base 
pressure  was  affected  by  Re. 

(2)  Effect  of  M:  The  experiment  done  by  Zarin^1^!  showed 
that  when  M£4.5,  base  pressure  decreased  as  M  increased. 

However,  when  M^4.5,  base  pressure  rose  with  increasing  M.  This 
shows  that  base  pressure  varies  according  to  different  patterns 
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at  hypersonic  and  supersonic  speeds. 

(3)  Effect  of  apex  angle  and  dullness:  As  the  dullness 
increases,  the  specific  base  pressure  p^/p^,  increases.  When  M 
and  dullness  remain  unchanged,  base  pressure  increases  with 
increasing  apex  angle. 

The  major  variables  affecting  base  pressure  are  introduced 
above.  Their  individual  effect  seems  to  follow  a  specific 
pattern.  However,  it  is  very  hard  to  analyze  the  combined 
effects  of  these  parameters.  Some  patterns  may  be  obtained  if 
shoulder  parameters  such  as  and  Re^  are  used  to  treat  the 
experimental  data  on  base  pressure  because  external  should  flow 
near  the  base  and  should  directly  affect  base  pressure.  A 
detailed  summary  was  made  on  this  subject  in  reference  [12] . 

Theoretical  analysis  on  base  flow  is  still  based  on  the 
model  proposed  by  Chapman  et  al.  The  Chapman-Karst  method  is 
semi-empirical  which  agrees  with  experimental  results  at  low 
speeds.  It  is  widely  used.  However,  its  disadvantages  are 
obvious  at  high  speeds.  The  boundary  layer  integral  method 
developed  later  by  Lees  and  Reeves t ^ ] , [14 ]  more  successful. 
It  does  not  rely  on  empirical  constants.  This  method  is  also 
different  from  the  conventional  integration  method.  Due  to  the 
fact  that  the  interaction  between  a  viscous  layer  (free  shear 
layer)  and  the  external  inviscid  flow  must  be  taken  into 
consideration  in  studying  base  flow  and  a  new  unknown  variable 
must  be  introduced  as  a  section  parameter  in  choosing  the 
velocity  section,  it  is  necessary  to  add  more  equations.  As  a 
result,  we  have  the  integral  equations  of  continuity,  momentum 


and  moment.  The  three  unknowns  are  the  external  Mach  number  Me, 
momentum  loss  thickness  6  and  section  parameter  a.  The  equations 
are: 

dMr  N,  d6  <V2  da  -V, 

dx  ~  D  ’  dx  ~  D  ’  dx  ~  D  (A) 

where  D,  N^.,  N2/  and  N3  are  functions  of  Me,  0,  a  and  x.  This 
set  of  differential  equations  does  not  seem  unusual.  However, 
when  we  integrate  over  the  direction  of  flow  from  an  initial 
cross-section,  it  is  not  possible  to  consider  the  effect  of 
downstream  motion  on  that  upstream.  In  reality,  at  D=0,  the 
integration  cannot  proceed.  It  is  necessary  to  repeatedly  modify 
the  initial  value  by  iteration  to  continue.  At  D=0,  it  is  the 
neck  of  the  wake,  which  is  a  critical  point.  The  critical  point 
theory  has  been  used  to  solve  base  flow  for  years.  However,  this 
method  is  limited  to  a  dull  body  with  relatively  low  M  at  the 
shoulder  L12J . 

The  work  done  by  Zhar  Wenlong  et  al  [15],  [16]  incxuded  in 
engineering  calculation  of  turbulent  flow.  The  base  heat  flow 
analysis  might  be  on  the  rough  side.  To  treat  this  problem 
better,  it  is  necessary  to  determine  the  mean  heat  flow  with  the 
aid  of  energy  conservation  after  the  relations  among  parameters 
in  the  re-circulation  zone  are  understood.  Meng's  straight  line 
method  [I7]  was  also  used  to  calculate  the  base  heating  problem. 


III.  Numerical  Method 


Many  authors 1 “  1^2]  USed  numerical  integration  of  the  N-S 
equations  to  study  base  flow  or  other  separation  flows. 

Numerical  methods  can  approximately  be  classified  into  three 
major  types:  1.  finite  difference,  2.  finite  element  and  3. 
integration  relation  or  straight  line  method^7!.  We  believe 
that  although  C.  Taylor  et  al  recently  used  finite  element  method 
to  calculate  an  incompressible  flow  field,  the  majority  of 
numerical  computations  of  base  flow  was  done  by  using  a  finite 
difference  method^2-*!'  l24^  .  In  the  following,  we  will  discuss 
incompressible  and  compressible  flows. 

(A)  Solving  N-S  equation  for  Incompressible  Flow 
The  N-S  equations  expressed  in  original  variables  (v,  p)  can 
be  written  as: 

5 -VC  =  -  -L  Krad  p+v  v,d  (la) 

div  v=0  (lb) 

Because  it  is  not  a  single  initial  value  problem  and  the 
time  term  9p/ 3t  is  not  included  in  (lb),  the  computation  is 
somewhat  difficult.  In  this  case,  the  continuity  equation  and 
equation  of  motion  are  not  simultaneous  equations  with  respect  to 
variables  v  and  p.  The  continuity  equation  imposes  constraints 
to  the  variable  v.  It  is  a  problem  to  be  addressed.  Earlier, 
the  constraint  div  v=0  was  relaxed.  In  addition,  D=div  v.  The 
equation  of  motion  becomes: 


The  difference  equations  of  (la)  and  2)  are 


m 

i  Km. 


a 

i 

s 


(3) 
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(4) 


where  L  is  the  difference  operator  for  all  terms  except  for  the 


time  and  dependent  terms  in  (la).  <$xx  and  <$yy  are  the  second 


kn+l 


order  difference  operator.  Now,  D^j  =  0  is  used  to  meet  the 
constraint  condition.  Equation  (4)  becomes  the  equation  for  p/p 
in  the  nth  layer,  i.e. 


<<*.. +  a„>  </>//» -  q:  . , +» «j„ + a.,)  z>:  , 


(5) 


The  computation  procedure  is  as  follows.  Based  on  the  known  ^ j 


and  D^j,  Rirj  is  determined  by  equation  (5).  Then,  (p/P)i,j  is 


-*n+l 


determined.  vi, j  is  then  calculated  by  substituting  the  results 
into  equation  (3).  At  this  time,  div  vn+^0.  The  final  result 
is  obtained  by  repeated  iterations.  The  equation  of  continuity 
is  met  by  repeated  iterations.  Later,  an  artificial  compression 
method  was  developed^] .  ^  time  term  is  introduced  into  the 

continuity  equation: 


dp  ,  .  - 

~dy  ~c’  div  v  -0 


(6) 
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The  above  equation,  although  there  is  no  physical  meaning,  is  /58 
consistent  with  the  equation  of  continuity  for  a  compressible 
flow  of  c  interpreted  as  the  sound  propagating  speed.  The 
value  of  c  must  be  chosen  to  ensure  that  the  solution  is 
convergent  to  become  a  steady  solution.  Recently,  another  method 
was  used  in  reference  [26].  v  and  w  were  used  as  two  variables. 

The  two  equations  are: 


da 

dt 


+  v  •  V  a 


V  ‘q  =  -  V  *  £ 


(A) 


In  order  to  make  v  and  w  satisfy  the  condition  that  their 
divergence-  is  zero,  a  "potential  correction"  is  added  to  each 
time  layer. 

(B)  Numerical  Solution  of  N-S  Equations  for  Compressible 

Flow 

As  early  as  1965,  Brailovskaya presented  a  two  step 
explicit  difference  scheme  with  second  order  accuracy  based  on 
unsteady  N-S  equations.  Allen  and  Cheng  Singyi^19!  improved  it. 
Viscous  terms  were  treated  by  using  the  method  similar  to  the 
Dufort-Fronkel  treatment  to  eliminate  the  time  limitation  due  to 
stability  requirement.  Later,  Brailovskaya ^ 27 1  took  the 
treatment  in  reference  [19]  into  consideration  and  introduced  an 
"artificial  viscosity  term"  with  second  order  accuracy  in  the 
equation  of  continuity.  It  has  a  "numerical  smoothing"  effect  on 
the  region  where  the  flow  field  gradient  varies  greatly.  In  the 
smooth  region,  this  term  is  near  zero. 

During  the  process  of  computing  a  steady  problem  based  on 


133 


-r. 


unsteady  equations,  the  results  obtained  at  any  instance  do  not 
reflect  the  actual  physical  process.  Only  when  t-*30,  it 
approaches  the  steady  solution.  The  advantage  is  that  the  amount 
of  computer  time  required  to  obtain  a  convergent  solution  can  be 
reduced.  Dorodnitsyn used  an  "artificial  unsteady"  method  to 


improve  the  computing  process,  in  view  of  the  concept  in 
reference  [28] ,  reference  [29]  introduced  an  adjusting  factor  in 
solving  the  N-S  equations  for  a  two-dimensional  base  flow.  The 
results  showed  that  the  rate  of  convergence  was  improved. 

The  MacCormack  scheme  is  famous  due  to  the  fact  that  it 
successfully  calculated  an  inviscid  shock  wave  containing  flow 
and  the  interference  between  the  shock  wave  and  the  boundary 
layer.  In  1969,  MacCormack' s  scheme^30!  was  an  explicit  two  step 
format.  The  step  length  is  seriously  limited  by  the  stability 
condition.  Later,  he  introduced  a  mixed  algorithm^31!  which 
included  the  MacCormack  explicit  scheme  and  the  corrected  time 
splitting  plan.  The  effectiveness  of  the  computation  was 
apparently  improved.  The  expanded  MacCormack  scheme  allows  the 
use  of  any  mesh  to  approximate  any  series  of  equations  in  a 
physical  plane. 

The  following  discussion  is  focused  on  the  plan  proposed  by 
MacCormack  in  1981^3^^.  This  is  a  conditionless  two  step  process 
starting  from  the  conservative  N-S  equations: 


c tu_ 

dt 


*F  ,  &G  =lF'  + 

dx  dy  dx 


±6. 

dy 


(7) 


difference  is  taken  in  two  steps: 
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where  I  A  I  and  |b  |  are  matrices  with  positive  eigenvalues  X  |Aj  and 
Xj  3  | •  They  are  related  to  Jacobi  matrix.  A=Sf/3U  and  B=3G/9U. 

In  addition,  X |A  |  and  X  |B|  include  viscosity  effect.  The  first 
terms  in  (8a)  and  (8b)  are  in  the  explicit  MacCormack  format.  The 
second  term  is  totally  implicit.  The  coefficient  matrix, 
however,  is  doubly  diagonal.  In  order  to  facilitate  direct 
computation,  the  conventional  three  diagonal  chasing  method  for 
the  implicit  scheme  is  avoided  to  save  the  computational  load  in 
every  step.  In  summary,  the  time  step  is  greatly  widened  by  (8a) 
and  (8b).  The  Courant  number  could  reach  102-10^.  However,  we 
must  point  out  two  things:  when  At  is  chosen  to  satisfy  the  CFL 
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condition,  the  scheme  will  automatically  be  transformed  to  the 
explicit  algorithm.  Second,  it  is  very  important  to  add  an 
approximate  dissipating  term  in  the  process.  When  it  is  used  to 
calculate  flow  in  the  base,  special  precautions  must  be  taken  in 
treating  the  boundary  condition. 

One  of  the  major  difficulties  encountered  in  the  implicit 
method  is  how  to  treat  the  non-linear  convection  term.  The  Beam- 
Warming  implicit  factor  decomposition  method^34!  as  expanded  by 
Steger^33^,  where  the  non-linear  convection  term  of  the  (n+l)th 
layer  is  calculated  based  on  the  linear  function  established  by 
the  Taylor  series  expansion  of  the  (n+l)t^1  layer  and  the  nfc^ 
layer,  can  be  used  to  calculate  the  viscous  flow  field  in  any 
shape.  Reference  [35]  showed  that  this  method  is  capable  of 
accelerating  the  stabilizing  process  in  calculating  the  base  flow 
as  compared  to  an  explicit  method.  A  great  deal  of  computer  time 
can  be  saved. 

Steger  pointed  out  that  the  non-linear  term  in  the 

conservative  Euler  equations  is  a  homogeneous  function  of  U. 

This  property  actually  allows  the  splitting  of  flow  vectors  F  and 
-► 

G  into  subvectors  by  similarity  transformation.  The 
corresponding  off-center  difference  direction  is  related  to  the 
eigenvalues  of  3f/3u  and  3G/3U.  The  difference  between  this  new 
algorithm  and  the  Beam-Warming  scheme  is  that  the  coefficient 
matrix  of  the  linear  equations  thus  derived  is  an  upper  (lower) 
triangle.  Its  inverse  can  be  directly  found.  This  means  that 
the  flow  field  has  a  subsonic  region.  It  also  provides 
reasonable  results.  In  reference  [36],  this  idea  was  used  to 
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solve  the  inviscid  portion  of  the  N-S  equations, 
part  was  computed  using  a  method  similar  to  that  developed  by 
Saul'yev^^] .  T^e  computation  is  as  follows: 


=  F-+6!F'+d ;  G'+dl  G')-^RHS 

[ 1 +<rrf>  +Trrki-M‘>] 


(9a) 
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where  A+  and  A“  are  the  matrices  formed  by  the  positive  and  /60 

negative  eigenvalues  of  A=3F/<tof  respectively.  B+  and  B“  are  the 

•+ 

matrices  formed  by  the  positive  and  negative  eigenvalues  of  B=9G/3 

it  it 

0,  respectively.  A  and  B  are  matrices  reflecting  the  effect  of 
viscous  terms.  The  asymptotic  off-center  difference  of  the 
convection  term  is  related  to  the  signs  of  the  eigenvalues  of  A 
and  B.  This  technique  can  improve  the  resolution  of  low 
viscosity  areas. 


Because  the  time  step  is  vigorously  limited  by  stability  in 


an  explicit  method,  even  when  various  steps  are  chosen  in 
different  directions  to  calculate  a  high  Re  number  flow,  the  cost 
is  very  high.  It  should  be  pointed  out  that  there  is  nothing  new 
in  solving  N-S  equations  with  an  implicit  or  hybrid  method. 
However,  it  is  rarely  used  to  solve  the  flow  in  the  base. 

The  hybrid  scheme  is  to  split  the  equations  into  several 
terms.  Each  term  is  approximated  by  an  explicit  or  implicit 
difference  operator,  with  respect  to  the  step  corresponding  to 
each  operator,  the  maximum  time  step  AtM  is  allowed  to  be  used. 

If  the  splitting  is  to  categorize  difference  operators  into  two 
kinds:  in  Class  A  AtM  is  equivalent  to  Ata  and  in  Class  B  At(1  is 
equivalent  to  Atb,  and  Ata  «Atb.  In  this  case,  the  most 
meaningful  thing  to  do  is  to  treat  Class  A  operators  as  implicit 
operators.  Due  to  the  fact  that  the  explicit  format  which  is 
rigorously  limited  by  the  stability  condition,  is  replaced  by  the 
implicit  approximation,  plus  the  appropriate  explicit 
approximation  remains,  it  is  much  simpler  than  the  full  implicit 
scheme.  It  is  particularly  favorable  for  solving  thin  viscous 
layer  problems.  This  method  can  be  further  simplified  because 
the  simpler  the  operator  becomes,  the  easier  the  implicit 
treatment  is.  In  reference  [31],  transversal  derivatives  were 
split  into  viscous-parabolic  and  inviscid-hypobolic  parts.  This 
treatment  improves  the  effectiveness  of  the  the  computation. 

In  conclusion,  numerical  methods  to  solve  N-S  equations  are 
under  development.  Their  effectiveness  and  optimization  process 
are  improving.  For  example,  Beam-Warming  combined  the  stable 
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linear  multi-step  method  with  an  approximate  factor  decomposition 
method  to  create  the  optimized  ADI  method.  This  type  of 
technique  is  growing.  The  goal  is  to  accelerate  the  rate  of 
convergence,  improve  accuracy  and  strengthen  effectiveness. 

IV.  Opinions  on  Developmental  Trend  of  Numerical  Methods  for 

Separated  Base  Flow 

(1)  Implicit  methods,  such  as  those  presented  by  Beam- 
Warming,  Briley-McDonald^38!  and  MacCormack,  are  unconditionally 
stable  when  solving  a  pure  initial  value  problem.  However,  it 
very  difficult  to  calculate  the  stability  with  boundary 
conditions.  Recently,  Beam-Warming t 39 1  cited  an  example  to 
explain  a  more  generalized  stability  analysis,  including  the 
effect  of  boundary  conditions.  Thus,  some  ambiguities  are 
explained.  The  theoretical  basis  for  this  type  of  stability 
analysis  is  the  theory  of  Kreiss-Osher-Gustaf f son.  The  results 
show  that  unconditional  stability  algorithm  agrees  with  the 
implicit  boundary  condition. 

(2)  In  order  to  meet  the  requirements  for  a  complicated 
flow  such  as  a  non-symmetric  profile,  especially  in  improving  the 
resolution  of  the  flow  field  near  the  base,  it  is  necessary  to 
rationally  arrange  the  nodes.  It  is  better  to  employ  a  body 
coordinate  system  or  a  multi-layer  treatment[40]-[43] . 

(3)  It  is  very  cumbersome  to  calculate  the  shock  wave  and 
secondary  shock  wave  for  a  complicated  supersonic  flow.  It  is 
very  difficult  to  determine  the  shock  wave  position  using  the 


capture  method.  The  best  approach  is  to  combine  these  two 

methods.  The  shock  wave  allocation  method  should  be  used  for 

strong  shock  waves  or  primary  shock  waves  of  influence.  The 

relatively  weaker  secondary  shock  waves  should  be  treated  with 

» 

the  wave  capture  method. 

(4)  We  should  pay  attention  to  combining  numerical  methods 
with  theoretical  analyses.  For  instance,  the  expansion  mechanism 
of  the  shoulder  should  be  properly  employed  in  the  base  flow 
calculation.  Three-dimensional  automatic  adjustment  of  nodes  can 
be  used  in  calculating  the  separated  flow  at  high  Re  numbers.  In 
a  turbulent  flow,  theoretical  analysis  and  experimental  data  can 
be  used  to  determine  the  step  length  of  nodal  points.  In 
summary,  it  is  very  meaningful  to  apply  the  results  obtained  from 
theoretical  analysis  to  improving  the  effectiveness  of  numerical 
methods. 
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Advances  in  Inverse  Boundary  Layer  Techniques  /63 

Li  Suxun 

(Beijing  Institute  of  Aerodynamics) 

At  the  present  moment#  the  complete  time-dependent  N-S 
equations  are  not  widely  in  a  wide  variety 

of  engineering  applications.  This  is  not  just  due  to  the  lack  of 
understanding  of  turbulent  flow#  lack  of  suitable  turbulent 
models#  or  limitations  imposed  by  economic  conditions  and  computing 
capacity  and  speed.  Instead#  the  physical  phenomenon  to  be 
reflected  in  reality  does  not  require  such  completeness# 
complexity  and  delicate  detail.  The  results  obtained  by 
simplified  N-S  equations  can  already  meet  the  requirements. 

Based  on  the  physical  problems,  it  is  possible  to  simplify  the 
equations,  boundary  conditions  and  initial  conditions. 

Consequently#  various  simplified  equations#  such  as  the  steady  NS 
equation#  low  dimension  number  NS  equation#  Euler  equation 
obtained  by  omitted  viscous  terms#  potential  equation#  linearized 
equation#  parabolic  NS  equations  obtained  by  neglecting  some 
viscous  terms  or  a  portion  of  the  equations,  and  boundary  layer 
equation#  can  be  obtained.  A  variety  of  numerical  methods  were 
developed  for  these  equations  to  form  a  new  branch  in 
computational  fluid  dynamics. 

One  of  the  most  frequently  encountered  typical  situation  is 
a  high  Reynolds  number  flow  field.  In  this  case#  the  viscosity 
effect  of  the  fluid  is  only  apparent  in  a  thin  layer  near  the 


surface.  The  scale  of  this  thin  layer  is  much  smaller  than  the 
characteristic  scale  of  the  flow  around  it.  Based  on  this 
physical  fact,  Prandtl  presented  the  boundary  layer  concept. 
Momentum  changes  normal  to  the  surface  are  not  considered. 

Higher  order  viscous  terms  in  the  NS  equations  are  completely 
omitted.  Consequently,  we  get  the  boundary  layer  equation.  The 
boundary  layer  equation  is  not  only  totally  different  from  the  NS 
equations  in  complexity,  its  mathematical  nature  is  also  altered. 
It  has  the  characteristics  of  a  parabolic  partial  differential 
equation.  When  the  flow  is  attached,  the  results  obtained  based 
on  boundary  layer  equations  agree  well  with  the  data  provided  by 
the  complete  NS  equations.  In  addition,  they  are  consistent  with 
experimental  results.  In  use,  the  requirements  for  computer 
speed  and  capacity  in  solving  boundary  layer  equations  are  less 
than  those  involved  in  solving  complete  NS  equations.  However, 
when  flow  separation  occurs,  singularity  appears  at  the 
separation  point  when  boundary  layer  equations  are  used. 
Goldstein^!  did  a  rigorous  study  on  the  properties  of  two 
dimensional,  incompressible  boundary  layer  equations  at  the  point 
of  separation.  He  pointed  out  that  boundary  layer  equations  are 
only  applicable  downstream  and  upstream  from  the  separation 
point.  At  the  separation  point,  the  equations  begin  to  exhibit 
singularity.  Furthermore,  the  specific  form  of  the  square  root 
singularity  was  derived.  Since  then,  Stewartson  (1958^1, 

1969  ^ ,  1970*4]),  Kaplun  ( 1967  l ) ,  Messiter  (1975I6J),  Cousteix 
and  Houdeville  (1981^)  also  did  further  work.  Dean  (1950^) 
proved  that  there  is  not  singularity  near  the  separation  point 


for  complete  NS  equations.  The  solution  is  still  an  analytical 
one.  Therefore,  the  singularity  is  introduced  by  simplifying  the 
NS  equations. 

The  advantages  of  boundary  layer  equations  are  well 
recognized  in  practice.  Is  it  possible  to  use  them  in  some 
separated  flow  field?  Because  flow  separation  is  frequently 
encountered,  it  is  desirable  to  search  for  a  new  technique  which 
not  only  preserve?  the  advantages  of  boundary  layer  equations  but 
also  avoids  the  difficulties  associated  with  separation  point 
singularity.  The  inverse  boundary  layer  technique  is  a  new 
method  developed  with  these  ideas  in  mind. 


Manuscript  received  on  June  22,  1984.  Revised  received  on 
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When  the  boundary  layer  concept  is  applied,  the  flow  field  / 64 
is  divided  into  an  inviscid  flow  region  and  boundary  layer 
region.  If  we  base  the  velocity  distribution  (or  pressure 
distribution)  at  the  outer  fringe  of  the  given  boundary  layer  and 
non-flow  conditions  at  the  surface  of  an  inviscid  flow  to  solve 
the  velocity  distribution  inside  the  layer  and  consequently 
derive  other  parameters  (such  as  surface  shear  stress, 
displacement  thickness,  momentum  thickness,  etc.),  then  the 
method  is  usually  called  a  forward  method.  As  described  before, 
at  the  point  of  separation,  the  equation  is  singular.  Thus,  the 
feasibility  of  progressing  downstream  is  limited.  Let  us  think 
if  the  boundary  conditions  are  changed  so  that  the  shear  stress 
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distribution  (or  boundary  layer  displacement  thickness 
distribution)  on  the  surface  is  given.  Furthermore,  its 
transition  near  the  separation  point  is  smooth  so  that  the  outer 
edge  boundary  layer  velocity  also  becomes  an  unknown.  It  is  then 
possible  for  boundary  layer  equations  to  exhibit  no  singularity 
at  the  point  of  separation.  This  is  called  the  inverse  method. 

It  not  only  preserves  the  special  features  of  boundary  layer 
equation  but  also  can  be  used  to  describe  the  characteristics  of 
a  flow  field  with  a  small  separation  zone.  Therefore,  it  has 
received  a  lot  of  attention.  Similar  to  the  forward  method,  the 
inverse  method  is  making  progress  in  two  directions: 
differentiation  and  integration. 

A  review  of  the  history  of  the  inverse  method  can  predict 
future  development  trends,  in  1966,  Catherall  and  Mangier ^ 
first  introduced  this  concept.  Moreover,  they  employed  the 
inverse  boundary  layer  integral  equation  technique  to 
successfully  pass  through  the  point  of  separation.  In  1968, 
Reyhner  and  Flugge-Lotz t10!  published  their  results  of  the  study 
on  the  interference  between  a  shock  wave  and  the  boundary  layer 
in  a  compressible  laminar  flow.  Differential  equations  were 
solved  with  the  inverse  method.  They  proposed  the  FLARE 
assumption  to  resolve  the  instability  problem  in  numerical 
computation  caused  by  the  reverse  flow  in  the  separation  zone. 

It  should  be  pointed  out  that  the  numerical  results  on  laminar 
flow  separation  bubbles  in  a  paper  published  by  Briley^^-^  in 
1971  by  using  steady  two-dimensional  incompressible  NS  equations 
quantitatively  described  the  flow  field  in  the  re-circulation 


zone.  The  inverse  in  boundary  layer  thickness  after  separation 
not  only  serves  as  the  basis  for  comparison  for  the  inverse 
technique  but  also  illustrates  physically  that  the  boundary  layer 
concept  is  still  applicable.  Keller  and  Cebeci(12l  in  1972 
compared  their  work  with  Briley's  results.  They  specified  the 
shear  stress  distribution  on  the  surface,  began  with  boundary 
layer  differential  equations,  and  used  the  FLARE  assumption  to 
obtain  successful  results.  In  addition,  Horton  (1974^^1), 
Klineberg  and  Steger  (1974  ^1)  gave  examples  to  calculate  two- 
dimensional  laminar  incompressible  separated  flow  fields. 
Different  difference  schemes  were  given  for  attached  flow, 
separated  flow  and  mixed  region.  Carter  (1974  (^J),  and  Wornom 
(1975 also  presented  corresponding  difference  schemes  and 
obtained  results  for  laminated  flow  with  separation  and  re¬ 
attachment.  Carter  published  computational  results  on  the 
interference  between  a  viscous  flow  and  an  inviscid  flow  in 
1975(17],  By  1978,  a  unified  Crank-Nicolson  difference 
scheme(l®l  suitable  for  both  forward  and  inverse  methods  was 
introduced.  It  is  applicable  to  two-dimensional  incompressible 
flows  and  compressible  laminar  and  turbulent  flows.  Furthermore, 
in  1979-1981  this  method  was  used  in  a  transonic  flow  field  with 
an  interference  zone  between  the  shock  wave  and. the  boundary 
layer.  A  satisfactory  iteration  method  was  resulted^'2®} . 

Kuhn  and  Nielson  (1973  121] )  calculated  the  flow  field  in  the  rear 
part  of  an  axisymmetric  body  with  jet  exhaust.  Kawai  (1977 t22l) 
presented  examples  of  calculating  the  compression  angle  in  a 
compressible  laminar  flow.  Williams  (1975(2^),  Whitfield  and 


Swafford  (1975 124] ^  a^so  obtained  successful  results.  Cebeci 
(1975^5]  f  1976  (26] f  1979(2?])  transformed  the  boundary  layer 
equations  into  various  forms.  In  addition,  the  accuracy  was 
improved  by  iterations  along  the  direction  of  flow.  Satisfactory 
results  were  obtained  in  a  variety  of  small  separation  regions. 

He  not  only  introduced  an  algebraic  turbulent  model  but  also  used 
the  complicated  two-dimensional  turbulent  model  in  the 
computation.  It  was  successfully  used  in  internal  flow  and 
trailing  flow  calculations.  Until  1983,  Agarwal  et  al(2®l  still 
based  on  a  two-equation  turbulent  model  to  solve  the  two- 
dimensional  and  axisymmetric  separation  flow  in  a  transonic  flow 
by  the  inverse  technique. 

On  the  other  hand,  many  worthwhile  results  were  obtained  by 
using  this  inverse  technique  to  integrate  boundary  layer 
equations.  Klemp  and  Acrivos  (1972 (29])f  East,  Smith  and 
Merryman  ( 1977  1 30 ] )  calculated  flow  fields  with  small  separation 
zones  based  on  Green's  roll-in  method.  Whitfield  (1978 1311) 
improved  the  velocity  model  in  the  separation  zone  and  obtained 
results  on  compressible  turbulent  layer  integral  equations  by 
using  the  inverse  approach.  Swafford  ( 1979  1 32 1  >  introduced  the 
analytical  expression  for  the  velocity  model  inside  a  two-  /65 

dimensional  turbulent  separation  zone.  Kline  (1981 ^ 3 ^ 3 )  used 
Cole's  wall-wake  law,  their  separation  criteria,  and  experimental 
data  to  obtain  various  results.  Based  on  various  methods 
published,  we  can  see  that  it  is  easier  to  be  successful  in 
calculating  a  viscous  or  inviscid  flow  by  using  the  inverse 
method  for  integral  boundary  layer  equations. 


In  practice,  the  iterative  solution  of  the  interference 

between  a  viscous  flow  and  an  inviscid  flow  receives  a  lot  of 

attention.  Ainarante  and  Chang  (1979  conducted  a  calculation 

on  the  interference  near  the  separation  bubble  of  a  laminar  flow. 

Kwon  and  Pletcher  (1979 f35^),  Veldman  (1979 t36^,  1981 f 37 ^ )  showed 

the  coupled  solution  for  an  incompressible  laminar  flow  with 

strong  interference.  In  recent  years,  when  there  is  shock  wave  - 

boundary  layer  interference  in  a  transonic  flow  field,  the 

iterative  solution  of  viscous  and  inviscid  flows  has  been 

developed  rapidly.  The  interference  at  the  jet  exhaust 

calculated  by  Kuhn  (1979  ^®^,  1980^®^),  the  steady  transonic 

flow  field  around  the  wing  estimated  by  Collyer  (1979 t4®])  and 

the  iterative  calculation  of  three-dimensional  time-dependent 

*■ 

Euler  equations  and  turbulent  boundary  layer  equations  by 
Whitfield  (1980 f24l)  belong  to  the  category  of  iteration  based  on 
the  inverse  method  of  integral  boundary  layer  equations.  Similar 
to  earlier  methods,  it  is  required  to  provide  a  semi-empirical 
formula  for  frictional  drag  and  momentum  thickness  (or  shape 
factor).  Carter  ,  based  on  differential  equations,  used 

the  inverse  method  for  finite  difference  boundary  layer  equations 
coupled  with  iterations  involving  the  complete  potential  flow 
equations  to  obtain  results  on  the  separated  region.  The  method 
is  simple.  It  converges  rapidly  and  is  applicable  to  flow  fields 
with  strong  interference  regions.  In  the  iterative  process 
involving  a  viscous  and  an  inviscid  flow,  there  are  two  commonly 
used  methods.  1)  The  boundary  velocity  of  the  inviscid  flow  is 
allowed  to  be  perturbed.  2)  The  boundary  layer  position 


thickness  (or  its  equivalent)  is  allowed  to  be  perturbed.  There 
are  successful  attempts  either  way. 
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Our  understanding  on  three-dimensional  separated  flow  is 
getting  deeper.  Flow  visualization,  topological  analysis  and 
various  numerical  computations  directly  promoted  the  investiga¬ 
tions  on  the  three-dimensional  inverse  technique.  Work  by  Le 
Balleur  (1981 141J )  ,,  stock  (1980^42bf  Formery  (1982 143l )  r  Blaise 
(1982^44J),  Delery  and  Formery  ( 1983  1 45 ] )  reveals  some  early 
information  in  this  aspect. 

In  recent  years,  work  in. this  area  began  in  China^46'  47f 
48,  49].  It  limited  to  some  new  studies  on  two- 

dimensional  laminar  and  turbulent  flow. 

A  review  of  this  development  process  makes  us  feel  that  the 

* 

mechanisms  for  some  complicated  flow  motions  in  high  Reynolds 
number  flow  fields  are  not  yet  well  understood;  including  three- 
dimensional  separation,  influction  separation,  relation  between 
separation  and  vortex  formation,  and  turbulence  models  in  and  out 
of  separation  regions.  For  this  reason,  more  assumptions  are 
used  in  establishing  a  mathematical  model.  Despite  the  use  of 
various  equations  and  numerical  computation,  it  is  often 
difficult  to  obtain  satisfactory  results.  Therefore,  in 
exploring  the  flow  mechanism,  we  must  be  able  to  identify  various 
major  influencing  factors.  We  must  offer  some  methods  to  verify 
theoretical  models.  In  addition,  some  meaningful  results  can  be 
obtained.  It  is  required  that  these  methods  are  not  too 
complicated  for  engineering  applications.  Based  on  the  fact  that 
the  inverse  method  discussed  above  is  one  of  those  methods,  its 


/so 


future  is  very  bright  throughout  the  world. 
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Abstract 

The  inverse  boundary  layer  techniques  have  been  developed,  which 
can  be  used  to  calculate  the  boundary  layer  with  separation  in  subsonic, 
transonic  and  supersonic  flow  field. 

A  review  and  summary  of  inverse  boundary  layer  methods  are  present¬ 
ed  in  this  paper,  including  both  laminar  and  turbulent,  in  2D  and  3D 
conditions,  as  well  as  viscous-inviscid  interaction.  It  has  been  shown  that 
the  inverse  techniques  provided  simple  and  efficient  tools  for  the  predic¬ 
tion  of  the  boundary  layer  flow  with  separated  region  and  will  play  an 
jmportant  role  in  future  applications. 
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Abstract 

This  paper  is  divided  into  two  parts.  The  first  is  about 
experimental  observations  of  vortex  breakdown  in  three  aspects: 
vortex  breakdown  types,  breakdown  region  structure  and 
observation  of  leading  edge  separation  vortex  breakdown.  All 
kinds  of  vortex  breakdown  phenomena  of  delta  wing  separation 
vortices  visualized  in  the  tube,  water  tunnel  and  low  speed  wind 
tunnel  since  its  first  discovery  in  1957  are  included.  The 
difference  between  vortex  and  separated  vortex  breakdown 
structures.  We  point  out  more  complicated  factors  affecting  the 
breakdown  of  wing  separated  vortex.  The  second  part  is  about  the 
theoretical  study  of  vortex  breakdown.  The  current  breakdown 
theories  are  divided  into  three  branches;  i.e.,  finite  transition 
theory,  hydrodynamic  instability  theory  and  vortex  core  theory. 

The  development  processes  of  these  theories  are  summarized. 

Their  advantages  and  disadvantages  are  compared. 

Some  interesting  problems  are  presented.  The  results  of  our 
recent  calculation  on  the  separated  vortex  breakdown  of  a  slender 
wing  is  presented. 


I.  Introduction 


The  vortex  breakdown  phenomenon  was  first  found  on  the 
leeward  surface  of  a  delta  wing^J.  Lambourne^  successfully 
obtained  two  types  of  vortex  breakdown  patterns  on  the  same 
photograph  of  the  leeward  surface  of  a  delta  wing.  Because  this 
phenomenon  is  closely  correlated  to  the  aerodynamic  performance 
of  an  aircraft  at  large  attack  angles^],  it  immediately 
attracted  the  attention  of  experimentalists  and  theoreticians. 
Investigations  are  centered  around  vortex  breakdown  mechanisms, 
forecast  of  conditions  to  cause  breakdown  and  various  types  of 
breakdown. 

Most  of  the  experimental  studies  concerning  breakdown 
mechanisms  are  performed  in  low  speed  circular  wind  tunnels.  A 
vortex  generator-fan  is  used  to  introduce  the  vortex  into  the 
pipe.  As  compared  to  separated  vortex  generated  by  a  wing,  it  is 
even  easier  to  control  the  influencing  parameters  with  chis 
method.  Harvey ^  first  used  this  type  of  apparatus  co  locate 
bubble  breakdown.  A  series  of  experiments  afterwards  such  as 
those  reported  in  references  [5]  and  [6]  further  proved  that 
Reynolds  number  and  vortex  strength  are  the  two  most  important 
initial  parameters  controlling  the  breakdown  mode.  Different 
combinations  can  result  in  six  vortex  disturbance  modes.  The 
impact  is  that  two  breakdown  modes  can  be  created  on  the  wing: 
symmetric  mode  and  spiral  mode.  They  can  also  be  created  in  a 
circular  wind  tunnel.  This  proves  that  the  boundary  condition  of 
the  separated  vortex  is  primarily  dependent  on  axisymmetry. 


Vortex  breakdown  theories  can  be  divided  into  three  types: 
finite  transition  theory,  hydrodynamic  instability  theory  and 
vortex  core  theory.  Although  they  are  able  to  resolve  a  series 
of  individual  problems  associated  with  vortex  breakdown  and  in 
some  cases  arrive  at  the  same  conclusions,  there  is  still  a  lot 
of  controversy  about  them. 

Recently,  the  detailed  structure  of  the  breakdown  region  was 
measured  by  laser  technology  in  references  [7]  and  [8].  For 
instance,  reference  [7]  reported  the  bubble  breakdown  of  a  double 
circulation  zone.  If  it  is  real,  then  even  original  axisymmetric 
N-S  equations  solved  by  numerical  method  such  as  the  work  done  in 
references  [9] ,  [10]  and  [11]  cannot  truly  reflect  the  breakdown 
mode.  Defects  in  existing  mathematical  models  are  fully  exposed. 
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II.  Experimental  Observations  / 69 

1.  Vortex  Breakdown  Types 

In  Sarpkaya's  experiment,  besides  seeing  bubble  breakdown 
and  spiral  breakdown,  he  also  discovered  a  double  spiral 
breakdown  mode.  Faler  et  al^l  conducted  experiments  over  a 
wider  range  of  Reynolds  number.  On  the  low  Reynolds  number  side, 
they  found  the  remaining  three  major  disturbance  modes.  By 
gradually  increasing  Reynolds  number  or  vortex  strength,  they 
sequentially  divide  the  major  disturbance  modes  into  six  types. 
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They  are  labelled  as  "6"  (off-center  shear  of  vortex). 


n  5  n 


(double  spiral),  "4",  "3"  (elliptical  bubble  in  transverse  cross- 
section),  "2"  (spiral),  "1",  "0"  (single  tail  bubble,  double  tail 


bubble) . 


These  experiments  confirmed  that  the  major  disturbance  mode 
is  only  related  to  the  initial  Reynolds  number  and  vortex 
strength.  The  pressure  gradient  of  the  external  flow  only 
affects  the  position  of  the  breakdown.  Figure  1  shows  the  three 


breakdown  boundaries  obtained  by  Sarpkaya.  In  specific  Re  and  Rc 
range,  a  dual  steady  state  exists.  It  was  called  the  lag  zone  by 
Sarpkaya.  In  thi-s  region,  there  is  a  tendency  to  spontaneous 
transformation  from  one  mode  to  the  other. 

2.  Detail  Structure  in  Breakdpwn  Region 


m 


Before  laser  speed  measuring  technique  was  available,  the 


internal  structure  of  the  breakdown  region  could  only  be 


displayed  through  flow  field  visualization  and  theoretical 


analysis.  In  the  past,  it  was  believed  that  a  bubble  is  a 


circulating  region  with  two  stationary  points.  However,  based  on 


the  measurement  made  by  Faler,  there  are  two  circulating  regions 


inside  a  bubble  type  of  breakdown.  There  are  four  stationary 


points  on  the  center  axis  as  shown  in  Figure  2.  Between  the 


first  and  second,  as  well  as  between  the  third  and  fourth 


stationary  points,  there  are  two  circulation  zones.  The  entire 


area  is  in  a  pulsing  state.  Figure  2  shows  the  results  after 


taking  an  average.  Measurements  reported  in  reference  [8] 


suggest  that  a  spiral  breakdown  field  vibrates  periodically.  Its 
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principal  period  is  equal  to  the  precession  period  of  the  spiral 
vortex . 
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All  measurements  mentioned  above,  indicate  that  every  point 
in  the  breakdown  region  of  the  flow  field  is  obviously  unsteady. 
Furthermore,  the  entire  breakdown  position  also  drifts  with  the 
equilibrium  position. 


Figure  1.  Boundaries  of  Various  Vortex  Breakdown  Types 
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Figure  2.  Double  Circulation  Zone  in  Bubble  Breakdown 


3.  Observation  of  Breakdown  of  Leading  Edge  Separated  Vortex 

Wing  separated  vortex  first  breaks  down  downstream  from 
the  trailing  edget1^].  As  the  attack  angle  increases,  the 
breakdown  position  is  gradually  shifted  toward  above  the  wing. 

The  aerodynamic  performance  deteriorates.  Because  the  Reynolds 
number  is  pretty  high,  only  bubble  breakdown  and  spiral  breakdown 
can  be  observed ^  •  f13 1 .  The  breakdown  position  is  dependent  on 
the  combination  of  attack  angle  and  sweepback  angle.  Figure  3 
shows  the  results  obtained  in  a  water  tunnel Increasing 
sweepback  angle  delays  the  breakdown  and  increasing  attack  angle 
moves  the  breakdown  forward. 

The  two  breakdown  modes  observed  on  separated  vortex  are  / 70 
similar  to  those  found  on  vortex  generated  in  the  pipe.  The 
difference  is  that  in  the  former  case  the  bubble  type  does  not 
have  a  secondary  breakdown  and  the  precession  direction  of  the 
spiral  is  opposite  to  the  spin  direction  while  in  the  latter  case 
they  move  in  the  same  direction.  It  is  worthwhile  to  point  out 
that  in  the  latter  case,  Reynolds  number  has  an  obvious  effect  on 
the  breakdown  type  and  position.  In  the  former  case,  it  is  not 
as  apparent. 
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Figure  3.  Effect  of  Attack  Angle  and  Sweepback  Angle  on  Delta 

Wing  Separated  Vortex  Breakdown 

1.  experimental  theory 

Reference  [15]  employed  a  hot  wire  meter  to  measure  the 
velocity  distribution  before  and  after  the  breakdown  of  wing 
leading  edge  separated  vortex.  The  patterns  obtained  are  very 
similar  to  those  obtained  in  the  pipe.  The  distribution  is 
extremely  close  to  axisymmetric.  The  phenomena  leading  to 
separated  vortex  breakdown  is  much  more  complex  than  those  for  a 
circular  pipe  vortex.  One  of  the  problems  is  vortex  formation. 
The  initial  condition  of  a  pipe  vortex  can  be  simulated  by  a 
correlation  equation.  However ,  the  mechanism  for  the  formation 
of  separated  vortex  near  the  wing  apex  is  still  not  understood. 
Next,  a  separated  vortex  is  under  the  influence  of  a  complicated 
external  flow  field.  At  large  attack  angles,  it  is  very 
sensitive  to  the  skid  angle^16^.  In  supersonic  range,  separated 
vortex  will  not  break.  In  transonic  range,  the  initial  breakdown 


position  begins  to  jump  around 


III.  Theoretical  Studies  on  Vortex 

To  date/  investigations  on  vortex  motion  and  its  breakdown 
are  based  on  incompressible  N-S  equations: 

div  ipy )  =  0  (1) 


+  (V  V  )  V  —  — ^-grad  P  +  -^~  div  r 


(2) 


Various  simplified  models  are  derived  based  on  these  equations. 


1.  Wave  Propagation  Theory.  In  an  incompressible  inviscid 
axisymmetric  vortex  motion,  the  finite  transition  theoxy  applies. 
The  force  involved  is  the  restoration  force  of  a  spring.  Hence, 
a  transverse  wave  of  infinitesimal  amplitude  will  propagate  along 
the  axis  against  the  flow  under  a  small  perturbation. 

Under  the  steady,  incompressible,  inviscid  and  axisymmetric 
assumptions,  equations  (1)  and  (2)  are  simplified  into  flow 
functions  in  cylindrical  coordinate  form: 


1>,.- r-V.  +  fc.-r*  H'  (y>)  -/'  (ip) 


(3) 


By  letting  y=r^/2  and  expanding  for  small  perturbation,  we 
get  a  second  order  differential  equation  for  ^(y) . 
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Together  with  the  boundary  condition  <4(0)  =<f>(a)  =o,  we  have  the 
Sturm-Liouville  system.  Based  on  Benjamin's  definition,  if  all 
eigenvalues  Ta  >  0,  then  the  flow  is  super  critical*  if  an  eigenvalue, 
such  as  vl  is  negative,  the  flow  is  subcritical.  y}«0is  the  critical  state. 

Squire^17!  first  found  the  condition  t*  *0.  He  believed  that 
the  infinitesimal  amplitude  wave  propagates  from  downstream  to 
upstream.  However,  when  it  reaches  the  critical  state,  a 
stationary  wave  of  infinite  wavelength  exists  in  the  vortex, 
blocking  the  propagation  of  the  perturbation  wave.  Hence,  the 
energy  accumulated  causes  the  breakdown  of  the  vortex. 

Nevertheless,  reference  [18]  proved  that  the  stationary  wave 
velocity  is  pointing  downstream.  Hence,  perturbation  energy  can 
only  spread  downstream.  Equation  (4)  has  two  mutually  conjugate 
solutions  corresponding  to  the  supercritical  and  subcritical 
states.  If  we  define  the  flow  force  (momentum  flow)  as: 

5  =  2wj.  r(p+  pu’)dr 

The  supercritical  flow  force  upstream  is  less  than  the 
subcritical  flow  force  downstream.  In  order  to  balance  the  flow 
forces,  breakdown  must  occur  in  between.  Moreover,  excess  energy 
is  consumed  by  the  broken  wake.  Hence,  vortex  breakdown  is  a 
finite  transition  from  a  supercritical  state  to  its  conjugate 
subcritical  state. 

Bossel^1^  overcame  the  defects  in  Benjamin's  model.  He 
assumed  that  the  flow  is  a  rigid  rotation  upstream.  Under  near 
critical  conditions,  a  series  solution  of  bubble  breakdown  was 
obtained.  Liebovich  et  al^°]  also  obtained  bubble  breakdown 


solution  by  employing  the  weak  non-linear  long  wave  theory.  He 
pointed  out  that  due  to  the  finite  transition  of  the  linear 
theory,  it  was  not  possible  to  obtain  critical  state  solution. 

It  is  a  non-linear  effect.  In  order  to  overcome  singularity  in 
the  critical  state,  the  flow  function  will  undergo  strange 
precession.  A  modified  Rorteweg-de  Vries  equation  was  obtained: 

A,  *=e(C,  AA.  +  C,  A..,)  +d'Ct(fA).  +  C>n  A 

It  is  an  initial  value  problem.  The  resulting  bubble  solution  is 
shown  in  Figure  4. 


Figure  4.  Single  Circulation  Zone  in  Bubble  Breakdown 
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Figure  5.  Vortex  Stability  Boundary 

1.  stable 

2.  unstable 

3.  stable 

Certain  additional  conditions  must  be  added  to  the  two 
methods  described  above  for  finding  the  near  critical  state  and 
critical  state  solution.  They  are  suited  for  studying  vortex  in 
a  circular  pipe,  including  numerical  solutions.  However,  they 
can  only  result  in  vortex  bubbles  with  a  single  circulation 
region. 

2.  Hydrodynamic  Instability  Theory.  Ludwieg^2^,  investigated 
the  hydrodynamic  stability  of  a  fluid  flowing  between  two 
concentric  circular  columns  and  then  expanded  his  results  to 
vortex  core  motion. 
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Let  us  consider  an  inviscid  flow  with  an  infinite  Reynolds 
number.  Furthermore,  let  us  assume  it  is  conical  and  the  vortex 
core  is  infinitely  slender.  Let  us  use  the  perturbation  method 
to  solve  simplified  equations.  The  perturbation  term  is  a  spiral 
perturbation.  Based  on  this,  we  get  the  stability  criterion: 
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Figure  5  shows  the  stability  boundary.  When  K<1.16  ,  vortex 
breakdown  occurs.  K  is  called  a  velocity  parameter.  This  theory 
believes  that  the  initial  spiral  perturbation  will  eventually 
lead  to  vortex  breakdown.  This  was  proven  in  reference  [22] .  It 
stresses  the  effect  of  non-symmetry  of  the  vortex  and  can  only  be 
used  to  analyze  spiral  breakdowns.  However,  this  theory  cannot 
take  the  axial  gradient  of  the  external  flow  field  into  account. 
On  the  other  hand,  spiral  perturbation  does  not  necessarily  exist 
in  the  upstream  flow  causing  vortex  breakdown.  Reference  [23] 
pointed  out  that  some  theoretical  analyses  even  proved  that  all 
incoming  flow  causing  breakdown  remained  stable  with  respect  to 
any  perturbation. 

3.  Quasi-cylindrical  Vortex  Core  Model.  Hall^24^  first 
established  the  quasi-cylindrical  model.  Based  on  analogy  in 
two-dimensional  boundary  layer  separation,  as  well  as  through  a 
series  of  numerical  calculations,  vortex  breakdown  occurs  when 
the  quasi-cylindrical  vortex  core  model  becomes  invalid. 

Inside  the  vortex  core,  v=eu  where  u  and  w  are  of  the  same 
order  of  magnitude.  If  we  further  assume  Re“^=0(c  ) ,  then  the 


vortex  core  model  equations  can  be  derived  from  equations  (1)  and 

(2)I25J: 
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It  is  parabolic  in  nature.  Hall  solved  these  equations  by  a 
difference  scheme.  References  (25).  [26]  and  [27]  obtained 
solutions  by  integration.  Mager^G]  employed  a  relatively  simple 
velocity  distribution  to  obtain  an  analytical  solution.  The 
axial  momentum  equation  (8)  and  radial  momentum  equation  (6)  are 
integrated,  respectively,  to  obtain 


'  a=(R  9,  —  K  ip)  ,'K  v 
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where  t =R(© j-0. 25  InR) .  Equation  (9)  is  a  parabolic  equation. 
When  =  0.163989,  the  upper  branch  is  connected  to  the 


lower  branch.  Equation  (10)  is  a  linear  equation.  0j*  can  be 
considered  as  the  critical  dynamic  loss.  The  intersect  of  curves 
described  by  equations  (9)  and  (10)  is  the  solution.  When  6^> 
©1*,  there  are  two  possible  solutions.  The  upper  branch  solution 
corresponds  to  the  supercritical  state  and  the  lower  branch 
solution  corresponds  to  the  subcritical  state.  The  "break¬ 
through"  from  the  upper  branch  to  the  lower  branch  corresponds  to 
Benjamin's  finite  transition;  i.e.,  bubble  breakdown.  The 
intersect  a  is  a  singular  point,  corresponding  to  spiral 
breakdown. 

References  [25]  and  [27]  relatively  systematically 
investigated  the  effect  of  Reynolds  number,  initial  vortex  core 
speed  and  circulation,  and  external  flow  speed,  pressure  and 
circulation  gradient.  Consequently,  the  breakdown  mechanism  was 
further  understood.  The  analysis  done  by  either  Benjamin  or 
Mager  on  the  breakdown  addressed  uniform  external  flow 
conditions.  As  for  a  non-uniform  external  flow  vortex,  even  if 
it  is  in  a  supercritical  state  upstream,  it  will  not  break  down 
if  the  external  flow  gradient  is  approximately  controlled. 
Reference  [28]  calculated  the  delta  wing  separation  vortex 
positions  at  a=65°,  70°  and  75°  by  using  a  double  layer  model 
which  agreed  well  with  experimental  results  as  shown  in  Figure  3. 
Recently,  some  authors  attempted  to  more  precisely  plot  the 
detailed  structure  of  the  breakdown  zone  by  taking  non-symmetry 
or  unsteadiness^®'  29] 


into  account. 


IV.  Conclusions 


Tube  vortex  and  leading  edge  separated  vortex  can  produce 
similar  axisymmetric  and  spiral  breakdowns.  Based  on  this  fact, 
axisymmetry  is  the  primary  factor  affecting  the  breakdown  / 73 

mechanism.  The  supercritical  and  subcritical  criterion  has 
already  been  widely  applied  to  vortices  inside  a  circular  pipe  to 
determine  whether  a  breakdown  will  occur.  Analytical  method  or 
numerical  solution  can  approximately  show  the  internal  structure 
of  a  broken  down  bubble.  However,  it  is  unable  to  calculate  the 
internal  structure  of  a  spiral  breakdown.  The  vortex  core  model 
accurately  predicts  various  breakdown  positions.  Furthermore,  it 
can  be  used  to  calculate  wing  separated  vortex  breakdown. 

Although  the  primary  property  of  the  vortex  is  axisymmetric,  it 
is  necessary  to  perfect  existing  mathematical  models  to  include 
non-ax i symmetry,  unsteadiness  and  compressibility  effects. 
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Numerical  Analysis  of  a  3-D  Separated  Flow  / 

Ma  Yanwen 

(Beijing  Institute  of  Aerodynamics) 

1.  Introduction 

A  missile  encounters  various  separated  flows  in  its  flight. 
An  important  technique  to  study  separated  flow  is  to  solve 
Navier-Stokes  equations  by  numerical  methods.  A  great  deal  of 
progress  has  been  made  in  the  numerical  simulation  of  two- 
dimensional  problems.  However,  there  are  difficulties  in  three- 
dimensional  flows.  In  addition  to  flow  complexity  itself,  we  are 
also  limited  by  computer  capacity  and  speed.  Computational 
methods  are  still  yet  to  be  improved. 

A  one-step  difference  scheme  was  used  in  reference  [1]  to 
calculate  a  simple  three-dimensional  separated  flow.  Their 
results  are  analyzed  here.  In  this  paper,  the  physical  model  is 
described  first.  The  calculated  results  are  then  analyzed.  In 
addition,  the  characteristics  of  the  three-dimensional  separated 
flow  are  discussed. 

2.  Physical  Model 

As  shown  in  Figure  1,  a  wedge  is  standing  vertically  on  a 
plate.  Over  the  plate,  there  is  a  uniform  supersonic  flow.  Near 
the  surface,  there  is  a  laminar  attached  layer.  The  wedge 
changes  the  direction  of  the  flow  (the  plane  of  symmetry  of  the 
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wedge  is  2=0) .  The  shock  wave  from  the  leading  edge  of  the  wedge 
interferes  with  the  attached  layer  on  the  plate.  In  this  work, 
the  interference  away  from  the  wedge  surface  is  analyzed  by  using 
numerical  results  (see  Figure  1). 


Figure  1.  .  Computation  Model 

1.  wedge 

2.  side  wall 

3.  shock  wave 

4.  bottom  wall 


Figure  2.  Upstream  Affected  Region  in  Interference  Zone 

1.  side  wall 

2.  shock  wave 

3.  upstream  affected  area 


The  N-S  equations  were  simplified  by  taking  into  account 
that  the  upward  velocity  component  perpendicular  to  the  bottom  is 
small  and  the  main  flow  is  in  the  x-direction.  In  order  to 
describe  the  flow  characteristics  in  the  interference  in  more 
detail,  the  y-direction  coordinate  was  transformed.  The 
simplified  equations  are  in  the  one-step  approximation  scheme. 

Manuscript  received  on  June  11,  1984,  revised  manuscript  received 
on  September  3. 

This  one-step  scheme  has  the  advantages  of  appropriate  accuracy,  / 
ease  of  finding  solution  and  good  stability.  In  order  to  reduce 
the  jumping  effect  in  the  numerical  solution,  a  smoothing 
parameter  is  introduced.  An  adjustment  factor  is  added  in  the 
difference  equations  to  accelerate  the  converging  process.  The 
boundary  conditions  given  are  shown  in  reference  [1], 

In  the  flow  field  calculated,  M  *  2.94  and  the  shock  wave 
angle  is  27.82°.  Two  sets  of  calculations  were  made:  one  for  Re6 
=  687. 5  and  the  other  for  Refi  =  3000.  The  mesh  points  in  three 
directions  are  25  x  31  x  31.  The  densifying  factor  in  the 
coordinate  transformation  is  b=25.  Thus,  there  are  approximately 
10  mesh  points  in  the  attached  layer.  A  stable  solution  can  be 
obtained  after  300  iterations. 


3.  Calculated  Results  and  Analysis 

Some  results  are  shown  in  Figures  2-11.  Figures  2-4  show 
the  first  set  of  results  and  Figures  5-11  are  the  second  set. 
Figure  2  shows  the  upstream  affected  area  due  to  interference 
between  the  shock  wave  and  the  attached  layer.  Figure  3  shows 
how  the  attached  layer  thickness  varies.  Based  on  this  we  know 
that  the  high  pressure  behind  the  shock  wave  makes  the  attached 
layer  thinner.  This  suggests  that  there  is  an  altitude  in  the 
interference  region  where  the  velocity  component  parallel  to  the 
inviscid  shock  wave  and  the  base  ug  is  smaller  in  front  of  the 
wave  than  that  trailing  the  wave.  This  is  one  of  the  reasons  why 
the  velocity  mode  changes  in  the  interference  region  as  we  will 
discuss  later.  Figure  4  shows  the  lateral  pressure  distribution 
at  1*20.  On  the  surface  j=lf  a  pressure  plateau  is  obvious. 
Figure  5  shows  the  variation  of  velocity  u  (velocity  component  in 
x-direction)  along  y-direction  on  the  cross-section  1=20  at  K=10, 
15  and  20.  K=10  is  located  behind  the  wave.  In  this  caser  u 

remains  in  a  Blasius  type  of  velocity  distribution.  K=15  is  in 
the  interference  region.  We  can  see  the  u  is  changing  mode. 

Figure  6  shows  y-direction  velocity  component  v  with  changes  in  y  at  three 
different  K  points  when  1=20.  We  can  see  that  in  the  wave  front  interference 
region  (K=20),  v  is  positive.  In  the  wave  tail  region  (K=10)  v  is 
negative.  According  to  the  pattern  of  u  and  v,  we  find  that  the  fluid  turns 
upward  near  the 

interference  region  and  then  deflects  downward  to  flow  through 
the  interference  zone.  The  z-direction  velocity  component  w  vs. 
y  curve  at  1=20  is  shown  in  Figure  7.  In  the  inviscid  region  in 
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front  of  the  wave  (K=20) ,  w=0.  This  means  the  flow  moves  along 
the  x-direction.  w  is  positive  in  the  viscous  layer.  The  fluid 
is  deflected  away  from  the  wedge  surface.  In  the  inviscid  region 
(K=10)  behind  the  wave,  the  fluid  moves  parallel  to  the  wedge 
surface.  The  flow  above  the  viscous  region  is  deflected  toward 
the  wedge.  In  the  viscous  region,  it  flows  away  from  the  wedge. 
This  effect  is  ever  more  apparent. 


\1  2.91.0 

Figure  3.  Variation  of  Attached  Layer  Thickness 

1.  shock  wave 

2.  Y-Z  plane  1=12 

3.  shock  wave 

4.  X-Y  plane 


Figure  4.  Pressure  Distribution  at  1=20 

1.  wedge 
. 2.  K-index 

Figure  8  shows  the  projection  of  velocity  in  the  normal  direction  / 77 
of  the  inviscid  shock  wave  un  against  y  at  R=10,  16,  and  22  when 
1=20.  At  K=16  and  22,  we  can  see  that  un  is  negative  near  the 
surface  •  This  indicates  that  the  fluid  flows 

backward  in  front  of  the  inviscid  shock  wave  above  the  viscous 
region.  On  the  bottom,  it  flows  forward.  The  un-v  pattern 
recorded  at  1=20  (see  Figure  9)  illustrates  the  formation  of  flow 
separation.  It  resembles  an  ellipse.  Its  long  axis  is  in  the 
circulation  zone  pointing  toward  the  surface.  The  high  pressure 
behind  the  shock  wave  makes  the  latter  half  of  the  circulation 
zone  stick  to  the  surface. 
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Figure  9.  un-v  Flow  Field  Pattern 

In  summary,  based  on  the  flow  pattern,  the  mass  in  the  upper 
half  of  the  viscous  layer  moves  downstream  along  the  x-direction  / 78 
and  is  slightly  deflected  upward  after  passing  through  the 
interference  zone.  Then,  it  is  deflected  downward  and  flows 
downstream  parallel  to  the  wedge  surface.  The  mass  in.  the  inner 
layer  moves  downstream  along  the  inviscid  shock  wave  plane 
direction  in  a  spiral  form.  Figure  10  shows  us  vs.  y  curves  at 
various  K  values  when  1=20.  ug  is  the  projection  of  velocity  in 
the  direction  parallel  to  the  inviscid  shock  wave  and  the  base. 

It  is  obvious  that  ug  has  modes  other  than  the  Blasius  mode.  It 
does  not  seem  to  be  caused  by  error  in  computation  and  be 
interpreted  physically.  This  change  is  caused  by  the  thinning  of 
the  attached  layer  due  to  the  high  pressure  of  the  latter  half  of 
the  interference  zone  and  the  spiral  motion  in  the  interference 
zone.  Because  ug  is  parallel  to  the  shock  wave,  the  value  of  ug 


is  equal  before  and  after  the  shock  wave.  In  addition,  because 

the  attached  layer  near  the  high  pressure  region  is  thin,  there 

at  MU  Qltftvde.  Ieve.1 

is  an  altitude^at  which  the  velocity  us  at  point  KA  is  less  than 
the  us  at  point  B  behind  the  inviscid  wave  front.  It  is  also 
known  that  the  mass  inside  the  interference  zone  moves  downstream 
in  a  spiral,  thus  a  low  velocity  ug  mass  point  at  A  will 
penetrate  the  shock  wave  front  and  flow  toward  point  B.  This 
will  slow  the  velocity  at  point  B  down  (with  respect  to  us  at  a 
point  far  away  behind  the  wave  at  the  same  altitude) .  This  is  a 
reason  why  the  model  of  ug  is  changed. 


Figure  10(a)  Velocity  Model  for  us  (1=20) 
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Figure  10(b)  Velocity  Model  for  uB  (1*20) 
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Figure  11.  Schematic  Diagram  for  Spiral  Motion 


1.  shock  wave 
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Exact  Vortex  Solutions  of  the  Navier-Stokes  Solutions  /80 

Wu  Jiezhi 

(Chinese  Institute  of  Aeronautics) 

Since  Rankine  presented  the  two-dimensional  vortex  structure 
with  vortex  core  in  1882,  for  over  a  hundred  years  only  a  few 
exact  vortex  solutions  of  NS  equations  were  found.  In  addition 
to  Oseen's  work  in  1912  which  extended  Rankine' s  vortex  to 
unsteady  viscosity,  the  next  progress  was  made  by  Burgers ^  who 
extended  the  work  to  three-dimensional  axisymmetric  vortex  and 
unsteady  vortex^2'2!,  and  by  Sullivan^  who  extended  the  work  to 
three-dimensional  vortex  with  double  cell  structure  (See  Figure 
1)  and  to  unsteady  cases .  Because  it  was  found  that  a  tornado 
has  a  'double  cell  structure,  Sullivan's  vortex  has  attracted  a 
lot  of  interest.  Nevertheless,  all  the  three-dimensional 
vortices  discussed  above  have  one  assumption  in  common';  i.e.,  the 
axial  velocity  w=2az.  It  becomes  infinity  when  z-~.  At  z=0, 
there  is  a  stationary  point  which  is  not  rational.  The  root  of 
the  problem  is  that  the  three  velocity  components  are  not 
properly  coupled  under  realistic  boundary  conditions.  As  a 
matter  of  fact,  the  presence  of  a  stationary  point  already 
suggests  that  a  solid  boundary  should  be  introduced.  With  regard 
to  the  interaction  between  a  semi-infinite  straight  vortex  and  an 
infinitive  plane,  because  of  the  lack  of  a  characteristic  length, 
the  solution  should  be  similar  to  the  conical  case.  Yih  et  al^ 
proved  that  only  the  conical  solution  can  provide  a  finite 
velocity  at  infinity.  Fence,  the  conical  solution  may  be  more 
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realistic  (except  at  the  origin)  than  the  axisymmetric  solution. 
It  is  a  natural  tendency  to  study  the  conical  vortex.  It  can  be 
considered  as  the  third  stage  of  studying  exact  vortex  solutions. 
In  this  paper,  the  author  vrll  briefly  introduce  the  studies  done 
in  this  area  (details  were  reported  in  reference  [7]). 


Figure  1.  Sullivan's  Double  Cell  Vortex 


It  was  experimentally  demonstrated that  only  the  vortex 
core  of  a  turbulent  vortex  has  a  tendency  to  expand  conically. 

It  is  not  true  for  a  laminar  vortex.  In  addition,  the  viscosity 
of  the  turbulent  eddy  near  the  core  is  approximately  constant. 

It  was  also  demonstrated  experimentally^9^  that  in  order  to  make 
the  vortex  radius  of  Sullivan's  solution  realistic  it  is  required 
to  choose  the  dynamic  viscosity  v  to  be  480  times  of  the 
molecular  viscosity.  This  suggests  that  we  should  consider  the 
turbulent  conical  solution  where  the  eddy  viscosity  is  variable. 
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Let  us  assume  that  the  flow  field  is  steady.  Based  on  the 


conical  axisymmetric  assumption,  the  velocity  components  in  the  a. 
and  e  direction  of  the  spherical  coordinate  system  (R,  a  ,e  )  can 
be  written  as 


(1) 


where  r =Rsina  is  the  distance  away  from  the  z-axis  (<*=0)  .  x=cosa 
is  the  only  independent  variable.  The  velocity  component  vR 
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can  be  expressed  by  a  continuous  function  F(x).  Let  us  assume  / 8 1 
that  the  eddy  viscosity  e  is  a  function  of  x  and  denote  e (x)  = 

1/ 2k  V(x) ,  The  constant  k  is  so  chosen  that 

Furthermore,  let  F  =  2v(.i-x,)f ,  Then,  after  substituting  equation 
(1)  into  the  spherical  incompressible  NS  equation,  we  can  derive 
the  following  integral-differential  equations: 


/>  i  p  _  K(x) 

J  J  v’Cl-x’)1 


n'  +  (2/+i^?>'+v-a^0-0 


(it) 


where 


K(x)  =  kiG(x)+  H{x), 

-tf2!(0) +  /’](*-  **)  +  Q(l-*)  -~H<  1) 
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Because  we  integrated  three  times  when  deriving  equation  (2a) , 
and  we  employed  a  no-sink  and  no-source  condition  at  the  vortex 
axis  x=l,  therefore,  there  are  two  integration  constants  P  and  Q. 
Let  us  further  specify  the  following  boundary  conditions. 

(1)  On  the  solid  wall  (x=0) ,  the  adherence  condition  is  met. 

(2)  On  the  vortex  axis,  the  regularity  condition  is  satisfied. 
These  two  conditions  lead  to 


u.y 
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n  =  f  =  f '  =  0  at  x=0 

A  =  0,  f  is  finite  at  x=l  (4) 

Based  on  equation  (2b) ,  fi=  0  is  always  a  solution  satisfying 
equation  (4) .  In  order  to  find  the  non-trivial  vortex  solution, 
we  must  introduce  (this  is  the  only  way)  a  normalization 
condition 

n ’  -  -1  at  x=l  (5) 

We  can  prove  that  n'(l)  =  -wR^  and  w  is  the  vortex  on  the  axis. 
Hence  equation  (5)  suggests  a  mechanism  which  maintains  the 
vortex  on  the  axis  invariant.  In  reality,  this  requires 
importing  energy  from  the  other  end  of  the  vortex  axis  to 
overcome  viscosity  loss.  This  is  equivalent  to  a  model  in  which 
the  interaction  between  the  tornado  and  the  ground  is  observed  by 
moving  the  tornado  forming  cloud  to  an  infinitely  high  altitude. 

Based  on  equations  (4)  and  (5) ,  we  get  Q=0  and  P=-M/k^  where 
M=H' (l)-H(l) .  In  addition,  we  can  derive  the  pressure 
distribution  formula.  Especially  at  x=0  and  1,  we  get 

at  z  =  0 


mj^  =  \2k-r: 
p  1  „ 


(6) 


at  z-axis 
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Up  to  this  point,  in  addition  to  the  parameter  k  which 
represents  the  Reynolds  number,  only  the  function  v(x)  remains  to 


be  defined. 

Before  we  discuss  the  complicated  situation  for  v,  let  us 
first  review  the  case  that  v=const.  Equation  (2)  was  first 
simplified  by  ( Golvdshtik  )  f into  the  following  elegant  form: 
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k‘G(x) 


v’a-*1)’ 
2/fl'«=0 


(7a) 

(7b) 


Now,  H(x)=0.  Therefore,  the  right  side  of  equation  (7a)  does  not 
include  f.  G(x)  is  still  expressed  by  equation  (3b).  Even  more 
significant  simplification  occurs  at  equation  (7b) ,  which  can  be 
directly  integrated.  However,  this  simplification  brings  another 
fundamental  difficulty.  Equation  (7b)  only  allows  n  to  vary 
monotonically.  Consequently,  either  0(0)  =0  or  0(1) =0  must  be 
abandoned.  Hence,  we  assume  0(1) =1  or  o(0)=l.  These  two  new 
assumptions  are  both  artificial.  Therefore,  either  P  or  Q 
remains  arbitrary  and  becomes  another  independent  parameter.  The 
assumption  n(l)=l  indicates  that  the  adherence  condition  remains 
valid  while  vortex  core  has  singularity.  Serrin^1^  established 
a  complete  mathematical  theory  for  this  case  and  proved  the 
presence  of  a  solution  with  respect  to  a  parameter  set  (k,  P) . 
With  respect  to  another  (k,P)  set,  the  solution  did  not  exist. 

It  was  found  that  when  the  solution  exists,  various  (k,P)  subset 
flow  fields  have  different  shapes.  There  are  single  cell 
solution  and  double  cell  solution.  On  the  other  hand,  the 
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assumption  Si(0)*l  means  that  the  regularity  of  the  vortex  core  is 
maintained  while  the  tangential  velocity  on  the  surface  is 
allowed  to  be  non-zero.  Yih  et  al^l  used  Serrin's  method  to 
study  this  situation  and  also  obtained  single  cell  and  double 
cell  solution  with  a  certain  (k,P)  set  (See  Figure  2) . 


Figure  2.  Laminar  Conical  Vortex  Solution  by  Yih 

Obviously,  although  these  solutions  are  simple  (numerical 
calculation  is  also  required  at  the  end) ,  they  are  not 
satisfactory  physically.  Furthermore,  the  significance  of  the 
parameter  P  is  not  clear.  In  reality,  their  results  proved  that 
if  steady  conical  vortex  does  exist  in  nature,  it  can  only  be 
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turbulent.  This  basic  difficulty  of  the  (  Golvdshtik-Serrin) 
equations  is  the  most  important  reason  why  we  investigated 
variable  viscosity.  It  is  obvious  from  equation  (2b)  that 
variable  vu )  is  a  necessary  condition  for  the  existence  of  a 
solution  of  equation  (4)  and  (5) . 

As  for  the  functional  form  of  v(x),  the  only  known 
limitation  is  the  symmetry  requirement  v'(l)=0.  Turner^2^  and 
Serrin^H]  believed  that  the  turbulence  in  the  tornado  has  a  / 83 

self-regulating  mechanism.  Hence,  the  problem  becomes:  What 
kind  of  form  should  *(*)  have  so  that  a  non-trivial  solution 
exists  for  equation  (2)- (5)? 

Let  us  first  assume  f  is  fixed.  Under  conditions  set  by 
equation  (4)  and  (5),  equation  (2b)  leads  to  a  singular 
eigenvalue  problem.  Through  some  expansion  of  existing  second 
order  eigenvalue  theory^13^,  we  proved  that  if 
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where  p  is  another  positive  parameter,  then  when  6(x;ji)$0  (0^x$l) 
there  is  only  a  trivial  solution.  When  6(x;^)  is  non-negative 
and  increases  monotonically ,  there  is  a  function  6  (x;/i)  which 
makes  [>an,>in+i)  have  a  finite  series  of  eigenvalues  x0,  x^,..., 
x i , . . , ,xn.  The  corresponding  eigen  function  n^(x)  has  i  zero 
points  between  (0,1). 

Next,  let  us  consider  the  entire  set  of  equations  (2).  By 
using  Schauder's  non-moving  point  theory  we  proved  that  with 
specific  parameters  (k,yu)  ,  there  is  a  set  of  solutions  (n,f) 
where  f>0  (0$x^l) .  Therefore,  it  is  a  single-cell  solution. 


These  existence  theories  are  sufficient  conditions.  Hence,  the 
possibility  of  a  double-cell  solution  is  not  excluded.  The  final 
check  will  require  numerical  calculation. 

It  can  also  be  proven  that  when  v(x)  is  monotonically 
increasing  between  0<x<e (e<<l)  and  then  remains  constant,  when  e- 
0,  our  solution  approaches  one  of  the  single-cell  solutions 
obtained  by  Yih*6^.  In  this  case,  the  simplified  equation  (2a) 
can  be  used  in  the  whole  region.  In  e<x<l,  equation  (2b)  is 
used.  There  is  a  boundary  layer  solution  for  n  in  0<x<e.  It  is 
only  related  to  6  (xjju)  and  independent  of  f.  This  boundary  layer 
solution  behaves  very  differently  from  that  of  Serrin's  solution 
at  x—0.  Therefore,  Serrin's  solution  is  not  meaningful  in 
maintaining  the  adherence  condition.  Yih's  solution  is  more 
practical. 

This  is  the  first  attempt  over  one  hundred  years  to  find  the 
exact  vortex  solution  of  the  NS  equations  (2) -(5)  which  satisfies 
all  the  physical  boundary  conditions  (except  the  origin)  by  using 
an  analytical  method.  This  investigation  using  a  simple  model 
fully  illustrates  the  complexity  of  the  motion  and  its  sensitivity 
to  the  boundary  conditions. 

Finally,  it  should  be  pointed  out  that  if  n(x)  has  i  zero 
points  between  (0,1),  then  there  are  i  co-axial  alternating 
vortex  loops  which  have  not  been  observed  in  practice.  They  may 


be  unstable. 
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Abstract 

The  characteristics  of  three-dimensional  separated  flow  are 
analyzed.  In  addition,  starting  from  the  fundamental  physical 
phenomenon,  based  on  the  results  of  three-dimensional  viscous 
separation  provided  by  the  numerical  solution  of  the  Navier- 
Stokes  equations,  as  well  as  the  experimental  data,  the  necessary 
conditions  for  three-dimensional  viscous  flow  separation  and  the 
definition  of  the  separation  streamline  on  the  surface  are 
presented. 

1.  The  Problem 

The  separation  of  attached  layer  is  an  important  phenomenon 
in  fluid  dynamics.  It  is  also  a  frequently  seen  phenomenon  in 
practice.  As  large  computers  are  being  developed,  a  great  deal 
of  progress  has  been  made  on  this  subject.  Earlier,  people 
thought  Tw=0  on  the  surface  is  the  separation  point  in  a  two- 
dimensional  viscous  flow.  Later,  Rott  employed  a  numerical 
method  to  solve  the  Navier-Stokes  equation  and  gave  the  flow 
field  of  a  plate  moving  in  a  static  fluid^.  His  results 
indicate  that  there  is  always  a  point  where  Tw=0.  However, 
there  is  no  separation  or  circulation  region.  In  fact,  it  is  not 


sufficient  to  consider  that  the  point  at  which  ^w=0  in  two- 
dimensional  viscous  flow  is  a  separation  point.  We  must  require 
that  a  separation  region  exists  behind  this  point  because  it  is 
possible  to  have  a  flow  where  tw>0  or  Tw =0  in  the  entire  flow^  . 
Moore,  Rott  and  Sear  proposed  the  famous  MRS  criterion:  in  a 
coordinate  system  moving  along  with  the  separation  point,  the 
conditions  for  flow  separation  in  a  two-dimensional  unsteady 
viscous  flow  are  Uq=0  and  (9u/3y)Q=0.  These  studies  resulted  in 
further  understanding  of  the  two-dimensional  separated  flow. 

However,  due  to  the  complexity  of  three-dimensional  viscous 
separated  flow,  many  phenomena  are  still  unclear.  In  recent 
years,  many  authors  proposed  various  criteria  for  the  separation 
of  the  three-dimensional  attached  layer.  Most  of  them  are  not 
precise.  In  order  to  make  the  following  discussion  easier,  four 
of  them  are  listed  below.  (1)  It  was  pointed  out  in  reference 
[3]  that  the  separation  streamline  of  the  three-dimensional 
attached  layer  is  the  boundary  of  the  circulation  zone  where  the 
fluid  cannot  penetrate.  (2)  References  [4,5]  pointed  out  that 
the  separation  streamline  is  the  envelop  of  the  limiting 
streamline.  Furthermore,  it  is  the  singular  streamline  of  the 
three-dimensional  attached  layer  equation.  (3)  Reference  [6] 
pointed  out  that  the  separation  streamline  is  the  limiting 
streamline  passing  through  the  point  Tw=0.  (4)  Reference  [7] 

pointed  out  that  the  necessary  conditions  for  the  three- 
dimensional  attached  layer  separation  are  rw.gradP=0  and  Uj.gradP> 
0,  where  P  is  pressure,  is  the  external  flow  speed  and  t*w  is 
the  surface  frictional  drag.  It  should  be  pointed  out  that  one 


of  our  interests  in  studying  the  conditions  for  the  separation  of 
viscous  flow  to  take  place  is  the  problem  itself.  On  the  other 
hand,  we  hope  to  determine  whether  separated  flow  exists  based  on 
the  flow  characteristics  near  the  separation  streamline.  It  is 
not  difficult  to  realize  that  the  first  three  criteria  cannot 
meet  this  goal.  Moreover,  they  are  somewhat  limited  and  are  not 
general  criteria  for  the  three-dimensional  separated  flow.  For 
example,  in  the  viscous  flow  around  a  sliding  cylinder  of  finite 
length,  there  is  attached  layer  separation.  However,  there  is  no 
envelop  of  the  limiting  streamline.  Nor  is  there  a  point  where 
'i^=0.  We  believe  the  latter  criterion  is  better.  Nevertheless, 
it  does  not  clearly  define  the  separation  streamline.  In  the 
following,  based  on  the  physical  properties  of  separated  flow, 
the  necessary  conditions  for  the  separation  of  three-dimensional 
attached  layer  and  the  definition  of  the  separation  streamline 
are  discussed. 

Manuscript  received  on  June  12,  1984.  Revision  received  on 
September  10. 

2.  Characteristics  of  Three-dimensional  Attached  Layer 
Separation  and  Its  Necessary  Conditions 

(1)  Definitions 

Let  us  consider  a  three-dimensional  viscous  flow  around  a 
curved  solid  wall.  Let  us  establish  a  perpendicular  coordinate 
system  (x,y,z)  where  z  is  the  normal  direction  of  the  curved 
surface,  x  and  y  are  on  the  tangential  plane  of  the  curved 
surface  (See  Figure  1) . 

i-jr 


In  order  to  facilitate  the  analysis,  let  us  introduce  the 
following  definitions: 

(i)  Wall  Frictional  Drag:  In  a  three-dimensional  viscous 
flow,  wall  frictioning  is  a  vector: 


(1) 


where  u  and  v  are  velocity  components  in  x  and  y  direction,  and 
juw  is  the  wall  viscosity  coefficient. 


Figure  1.  Schematic  Diagram  of  the  Coordinate  System 


(ii)  Limiting  Streamline:  We  realize  that  the  solid  wall 
has  no  flow.  There  is  no  streamline  on  the  wall  surface.  The 
streamline  immediately  next  to  the  wall  is  defined  as  the 
limiting  streamline.  The  slope  of  this  streamline  on  the  wall 
can  be  defined  by  the  following  equation: 


dx  u_ 
dy  v 


/f* 


(2) 


Let  us  expand  u  and  v  into  Taylor  series  at  z=0: 


u  = 


v  = 


• Z*-rO(z! ) 


•2:  +  0(2’) 


(3) 


By  substituting  equation  (3)  into  equation  (2)  and  taking 

equation  (1)  into  account,  the  following  expression  applies  to 

the  limiting  streamline  (z-»0): 

dx  ^t,,  (4) 

dy  r„. 

Based  on  equation  (4)  we  know  that  the  tangential  direction  of 
the  limiting  streamline  coincides  with  that  of  the  vector  *fw. 

(iii)  Accelerated  Flow  and  Retarded  Flow:  The  momentum 
equation  for  an  ideal  inviscid  flow  has  the  following  form: 

dW  (5) 

-jr~-2V  •®rad  P,p 

We  can  see  that  when  U.gradP<0,  the  velocity  is  in  an 
increasing  mode  and  when  U.gradP>0,  the  velocity  is  in  a 
decreasing  mode.  We  define  that  in  viscous  flow  that  U.gradP<0 
is  accelerated  flow  and  U.gradP>0  is  retarded  flow.  Here,  U  is 
the  velocity  of  the  point  mass,  P  is  the  pressure  and  p  is  the 
density. 

(2)  Properties  of  Two-dimensional  Attached  Layer  Separation 
and  Its  Criterion 

As  a  fluid  flows  around  a  solid  wall,  when  the  Reynolds 


number  is  sufficiently  high,  its  viscosity  effect  is  primarily 
exhibited  in  the  attached  layer  near  the  wall.  The  cause  of  flow 


separation  can  be  simply  described  as  follows:  due  to  frictional 
losses,  when  the  residual  kinetic  energy  of  the  fluid  in  the 
attached  layer  near  the  wall  is  not  sufficient  to  overcome  the  / 8 7 
pressure  to  let  the  fluid  flow  into  the  high  pressure  area 
downstream,  the  fluid  is  separated  from  the  wall.  At  the  same 
time,  the  mass  outside  the  attached  layer  still  flows  forward. 

Based  on  this  we  know  that  the  two  major  factors  causing  the 
separation  of  the  attached  layer  are  the  friction  drag  tw  and 
pressure  gradient  gradP  and  their  relationship  with  the  kinetic 
energy  of  the  fluid.  In  three-dimensional  flow,  when  the  fluid 
flows  along  an  arbitrary  curved  surface,  the  streamline  of  the 
ideal  flow  is  curved.  Moreover,  the  direction  of  the  pressure 
gradient  and  that  of  the  velocity  of  the  inviscid  ideal  flow  U-^ 
do  not  coincide  (See  Figure  2  in  which  the  x  and  y  direction 
velocity  modes  are  plotted) , 

Under  such  a  pressure  gradient,  the  fluid  mass  inside  the 
attached  layer  is  deflected  more  than  that  outside  the  attached 
layer.  On  top  of  it,  due  to  friction,  the  kinetic  energy  of  the 
mass  near  the  wall  is  gradually  decreasing  due  to  friction.  When 
the  kinetic  energy  of  the  mass  cannot  resist  the  forward  movement 
due  to  the  negative  pressure  gradient,  separation  takes  place. 

At  this  time,  the  projection  of  the  velocity  of  the  point  mass 
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near  the  wall  is  less  than  or  equal  to  zero  in  the  gradP 
direction;  i.e.,  U.gradPiO.  Hence,  the  necessary  condition  for 
attached  layer  separation  to  take  place  is  that  such  a  flow  area 
exists.  In  this  region,  the  flow  is  accelerated  inside  the 
attached  layer  and  it  is  retarded  outside  this  layer.  The 


m 


projection  of  the  boundary  of  this  region  on  the  wall  is  the 
separation  streamline.  Taking  into  account  that  the  tangential 
direction  of  the  limiting  streamline  coincides  with  tw,  the 
necessary  condition  for  the  separation  of  the  three-dimensional 
attached  layer  was  written  in  the  following  form  in  reference 
17]: 

"^.gradP=0,  U^.gradP  >  0  (6) 

Here,  U-^  is  the  velocity  vector  outside  the  attached  layer.  The 
separation  streamline  will  be  discussed  as  follows.  Because  grad 
P^O,  based  on  equation  (6)  we  can  see  that  there  are  two 
situations  where  the  first  condition  is  valid: 

(1)  |T,j  =  Oi.  +  T.j]'  ;=  0  ,  BP  i„-r„=  0  i  (2)  t  i. _L graaP 


Figure  2.  Schematic  Diagram  of  the  Velocity  Model 

We  believe  that  the  wall  streamline  is  the  limiting 
streamline  whose  tangential  direction  is  perpendicular  to  gradP. 


When  at  least  one  I  Twl  =0  point  exists,  the  separation  streamline 
is  an  approximation  of  the  limiting  streamline  passing  through 
this  point.  For  instance,  in  the  zero  attack  angle  viscous  flow 
around  a  slender  body,  the  pressure  gradient  goes  along  the 
surface.  Therefore,  the  separation  streamline  is  an  approxima¬ 
tion  of  the  limiting  streamline  perpendicular  to  the  object. 

Every  point  on  the  separation  line  corresponds  to  F£wl=0. 

However,  when  an  attack  angle  exists  for  the  same  flow,  two 
points  on  the  separation  line  correspond  to  |tw|=0;  i.e.,  the  two 
intersects  between  the  separation  streamline  and  the  plane  of 
symmetry.  When  there  is  nolfwl  =  0  point,  the  separation 
streamline  may  be  an  approximation  of  the  limiting  streamline 
perpendicular  to  gradP.  It  may  also  be  a  limiting  streamline. 

In  the  latter  case,  near  the  separation  streamline,  all  limiting 
streamlines  are  perpendicular  to  gradP. 
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Abstract 

Based  on  the  viscous/inviscid  interaction  model  and  the 
integral  boundary  layer  equations,  the  equations  for  two- 
dimensional  and  three-dimensional  incompressible  laminar 
viscous/inviscid  interaction  are  derived.  These  equations  can  be 
used  in  mild  cases  of  separation  and  re-attachment.  The  method 
and  examples  are  given  in  this  paper.  The  two-dimensional 
results  agree  with  other  computational  and  experimental  results. 
This  method  has  the  advantage  of  being  simple  and  fast'. 

1.  Introduction 

It  is  well  known  that  viscous/inviscid  interaction  in 
separated  flow  is  very  important.  Although,  in  principle,  this 
problem  can  be  solved  by  Navier-Stokes  equations,  however,  it  is 
a  very  expensive  and  time-consuming  way  in  practice.  In  recent 
years,  a  great  deal  of  progress  has  been  made  in  developing 
methods  to  calculate  viscous/inviscid  interaction  based  on 
boundary  layer  theory.  These  methods  include  direct  and  inverse 
techniques.  Goldstein  et  altl"3]  pointed  out  that  near  the  two- 


dimensional  separation  point  and  re-attachment  point,  direct 
boundary  layer  methods  will  have  singularity.  Thus,  direct 
methods  cannot  be  applied  to  separated  flow.  In  the  inverse 
method,  a  given  displacement  thickness  or  wall  shear  force 
distribution  can  be  used  to  avoid  singularity .  The  major 
characteristics  of  direct  and  inverse  methods  are  shown  in  Figure 
1. 
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Figure  1.  Schemes  of  Direct  and  Inverse  Methods  and  This 
Technique 

1.  inviscid  flow 

2.  boundary  layer  equation  (singularity,  if  strong 
interference  appears) 

3.  (a)  direct  method 

4.  initial  value 

5.  inverse  method  for  boundary  layer  equation 

6.  inviscid  flow 

7.  inverse  method 

8.  inviscid  flov* 

9.  interaction  pattern  and  boundary  layer  equation 

10.  this  method 


Stewartson[7]  introduced  the  three  layer  theory.  On  this 
basis,  Veldmanl8'S>l  presented  an  interacting  boundary  layer 
model.  In  this  model,  the  velocity  at  the  outer  fringe  of  the 
boundary  layer  Ue  is  not  known.  Instead,  it  is  a  part  of  the 
solution.  Furthermore,  he  assumed  that  Ue  is  a  linear 
combination  of  the  potential  flow  solution  Ue  and  the 
perturbation  velocity  Ue&*  due  to  viscosity,  i.e.: 

Manuscript  received  on  March  17,  1984.  Revised  manuscript 
received  on  September  5. 

ue  =  ue0  +  ue6*  M  /90 

He  employed  a  difference  method  to  calculate  the  two- 
dimensional  laminar  boundary  layer  with  separation.  Based  on 
Veldman's  interacting  model,  together  with  an  integral  form,  the 
authors!^0'!1!  presented  an  approximate  method  to  calculate  the 
two-dimensional  or  three-dimensional  incompressible  laminar  or 
turbulent  flow  with  separation  and  inviscid/viscous  interaction 
(see  Figure  1) .  This  method  will  be  described  in  detail  in  this 
paper.  New  examples  are  given.  For  ease  of  discussion,  only 
laminar  flow  is  involved.  Please  refer  to  references  [10,11]  for 


turbulent  flow. 
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2.  Computation  of  Two-dimensional  Interacting  Boundary 
Layer  With  Separation 


(1)  Integral  Equation 

The  two-dimensional  incompressible  laminar  boundary  layer 
integral  equation  can  be  written  as: 


dx  U.  dx  2 


where  CJe  is  the  velocity  at  the  outer  fringe  of  the  boundary 
layer,  0  is  the  momentum  thickness,  H  is  the  shape  factor,  and  Cj 
is  the  local  surface  friction  coefficient.  We  employed  Thwaites' 
method!^] #  The  parameter  X  is  defined  as: 


,  0*  dt/, 

V  dx 


where  v  is  the  dynamic  viscosity  coefficient.  In  addition,  the 
viscous/inviscid  interaction  pattern,  equation  (1) ,  is  adopted. 
Veldmanf®'^  pointed  out  that  in  many  cases  it  is  only  necessary 
to  use  the  approximate  expression  for  Ue  *.  If  we  adopt  the 
incompressible  thin  wing  theory, 


Here,  (xb,xe]  represents  the  important  area  for  interaction. 
Equations  (1)  -  (3)  form  the  basic  equations  for  the  two- 
dimensional  laminar  viscous/inviscid  interaction.  After 


mathematical  operations,  these  equations  can  be  written  as: 


dk  r~  dU,  r~  dd*  c 
-djr=F>i  di;=F 


(5) 


Here,  ,  F£,  and  F3  are  functions  of  boundary  layer 
parameters  (such  as  H,  6*,  Cj,  etc.).  The  details  are  shown  in 
reference  [10]  or  the  appendix  of  this  paper. 

(2)  Computing  Method 

The  initial  condition  of  the  above  equations  is  that  all 
parameters  of  the  initial  point  are  given.  Because  F^  contains  . 
the  unknown  8*,  iterations  are  required.  First,  an  assumed 
value  of  8*(x)  is  given  (usually  the  8*  distribution  of  the 
laminar  plate  boundary  layer  is  used) .  The  equations  can  be 
pushed  downstream  by  the  Runge-Kutta  method.  The  above  procedure 
is  repeated  until  8 *  converges.  The  iteration  equation  is 

+  a*'”  (6) 

Here,  6*  is  the  newest  value  of  6*  and  ©  is  a  relaxation 
parameter.  Usually,  super  relaxation  can  be  employed,  e.g., 

©=  1.5 , 


(3)  Calculated  Results  and  Discussion 

The  shape  calculated  -  a  "ditch"  is  shown  in  Figure  2.  The 
shape  of  the  ditch  is  given  by  the  following  equation: 

y “  —  0.03  sech  4CJf-2.o) 

leST 
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Figure  2.  Profile  Used  in  the  Calculation 

1.  equivalent  body 

2.  separation  point 

3.  re-attachment  point 

4.  interacting  region 

The  calculation  begins  at  x*l.  The  results  are  shown  in 

Figure  3.  The  results  of  6*  agree  very  well  with  Veldman's 

result^.  The  pressure  coefficient,  CD,  distribution  are 

Cf  * 

essentially  in  agreement  witty  the  results  obtained  by  Veldman^^ 
as  well  as  those  by  Carter  and  Wornomt13^ .  The  calculation  is 
rapid.  It  only  took  4  seconds  to  calculate  this  example  on  a 
VAX-11. 
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Figure  3.  Results  of  Laminar  Separation  at  the  Ditch 

1.  this  work 

2.  interaction 

3.  inviscid 

4.  inviscid,  this  work 

5.  interaction,  this  work 

6.  this  work 
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Figure  4.  Separation  of  Plate  Attachment  Layer  Under  Cylinder 


Another  example  is  to  calculate  the  plate  laminar  boundary 
layer  under  a  cylinder.  The  calculation  did  not  include  the 
effect  of  the  wake  behind  the  cylinder.  The  separation 
conditions  at  two  locations  are  shown.  Figure  4  shows  the 
separation  bubble  behind  the  cylinder  center.  Figure  5  shows 
that  separation  bubbles  can  appear  in  front  and  behind  the 
cylinder  center.  In  this  case,  the  cylinder  is  very  close  to  the 
plate.  The  example  demonstrated  the  ability  of  this  method  to 
deal  with  some  complex  flows. 


3.  Calculation  of  Three-dimensional  Interacting  Boundary 
Layer  With  Separation 

(1)  Integral  Method 

A  rectangular  coordinate  (x  y  z)  is  used.  In  this  work,  x 
is  the  free  flow  direction.  The  shape  of  the  object  is  given  by 
y=f(x,z).  The  momentum  integral  equations  of  the  three- 
dimensional  incompressible  boundary  layer  in  x  and  z  direction 
are: 


(7) 


Where  ®21r  ®22f  ^1  an{^  ^2  are  defined  in  the 

conventional  sense.  Ue  is  the  combined  velocity  at  the  outer 
fringe  of  the  boundary  layer.  Its  components  in  x  and  z 
direction  are  and  W^,  respectively.  Cf^  and  Cfz  are  the 
surface  friction  coefficients  in  x  and  z  direction,  respectively. 

Based  on  the  concept  of  Lighthill ,  the  three-dimensional 
equivalent  displacement  thickness  <3*  obeys  the  following 
equation: 
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The  viscous/inviscid  interaction  pattern  is  the  same  as  the 
two-dimensional  case;  i.e.,  equation  (1).  It  is  written  in  a 


component  form: 


Vi  =  vlo  +  vl6* 

Wi  =  Mlo  +  Wl6* 


By  using  the  thin  wing  theory  approximation 


Where  r  =  -^/(x-£)*  +  (z-rj)z.  £  and  r\  are  in  the  x  and  z 
direction,  respectively.  The  integration  region  covers  all  the 
viscous/inviscid  interaction  significant  region. 


(2)  Velocity  Cross-section  and  Equations 
In  order  to  solve  the  above  equations,  velocity  cross- 
section  and  auxiliary  equations  are  the  key.  Usually,  the 
velocity  cross-section  is  decomposed  into  flow  direction  and 
transverse  direction  components.  Coles pointed  out  that  the 
flow  direction  velocity  cross-section  of  the  three-dimensional 


flow  has  the  properties  of  the  two-dimensional  flow.  In  this 
work,  we  employed  Thwaites1  method  to  treat  the  flow  direction 
velocity  cross-section.  The  transverse  velocity  cross-section  is 
more  complicated.  There  is  no  general  expression  available  to 
date.  In  this  work,  we  used  Hager ' s 1  incompressible 
transverse  velocity  cross-section.  Thus,  all  transversal 
integral  thickness  is  related  to  the  momentum  thickness  in  the 
flow  direction  ^l*  Their  relationship  is  shown  in  reference 
[16]  . 

By  performing  a  geometric  transformation  from  the  streamline 
coordinate  system  to  the  rectangular  coordinate  system  and  using 
Mager's  expression  for  the  transverse  flow,  equations  (1),  (7), 

(8)  and  (9)  can  be  converted  into  the  following ^ 10^ : 


“  dx  ■  "  '  ' dx  U,  V  "-U.  "U,  ')  dx  ^U\ 

^  )&♦ 

¥jl  &A1-F  d9"  -6  F  dv+(a  F  w*  a  r 

u.  dX  dx  e“F',d7+V''F'.  ui~6“F'+fr;)~dt' 

+ 1±±.( f  v  +f  w  1  -  <: 

dKj-K  u  **l--dv>,  +  s 

dx  *lU'‘dx  _ar+5* 

™Li-K  U  dd*-dl1':  c 


Here,  a  is  the  angle  between  the  streamline  and  the  x-axis 
and  Y  is  the  angle  between  the  outer  fringe  velocity  of  the 
boundary  layer  and  the  frictional  force  line.  and  W-^  are  the 
velocity  components  at  the  outer  fringe  of  the  boundary  layer 
which  are  "floating"  as  unknowns.  They  are  determined  by  the 
viscous/inviscid  interaction.  V-,  and  W-,  correspond  to  the 
inviscid  velocity  components.  Ue  is  the  overall  inviscid 
velocity.  Because  the  right  side  of  the  equation  contains  the 
unknown  &*,  it  requires  iterations  to  solve  this  equation.  At 
each  mesh  point,  there  are  5  unknowns;  i.e.^i'  y,  V^,  and  6* 
There  are  five  equations.  The  boundary  condition  is  to  give  all 
parameters  at  the  upstream  Doundary  and  the  side  boundary  where 
the  gas  flow  enters  the  stream. 
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Figure  5.  Plate  Attached  Layer  Separation  Under  Cylinder  (II) 


Figure  6.  Plate  Attached  Layer  Separation  Under  Sphere 

(3)  Computational  Method 

The  procedures  to  solve  the  above  equations  are  as  follows: 
First,  assume  the  6*  distribution.  Usually,  a  laminar  plate  6* 
distribution  is  adopted.  Next,  calculate  the  partial  derivatives 
with  respect  to  z.  Here,  we  must  consider  the  direction  of  the 
transversal  flow.  On  the  side  boundary,  use  a  forward  or 
backward  difference  scheme.  Then,  use  the  Runge-Kutta  method  to 
proceed  in  the  x-direction.  Each  step  uses  the  new  value  of  6*. 
After  reaching  the  downstream  boundary,  it  is  re-calculated  until 
the  &*  distribution  becomes  convergent.  The  iteration  equation 
is  still  equation  (6) .  In  most  cases,  it  is  required  to  use  low 
relaxation,  such  as  £>  =  o.50 


(4)  Results  and  Discussion 

The  example  calculates  the  plate  boundary  layer  under  a 
sphere.  Similarly,  the  effect  of  the  wake  is  neglected.  Figure 
6  shows  the  local  friction  coefficient  distribution  along  the  x- 
direction.  We  can  see  that  the  computation  illustrates  the  local 
separation  case. 

As  mentioned  above,  &*  iteration  requires  low  relaxation 
which  is  different  from  the  super  relaxation  case  in  the  two- 
dimensional  case.  In  that  case,  only  one  velocity  component  must 
be  derived  from  &*.  In  the  three-dimensional  case,  however,  two 
velocity  components  must  be  formed  based  on  &*.  Hence,  these  two 
velocity  components  are  coupled.  Therefore,  in  the  three- 
dimensional  case,  the  number  of  iterations  required  far  exceeds 
that  in  the  two-dimensional  case. 

4.  Conclusions 

This  paper  introduces  an  approximate  method  to  calculate  the 
viscous/inviscid  interaction  with  separation  based  on  basic 
boundary  layer  assumptions.  The  results  indicate  that  the 
velocity  correction  Up  *  due  to  viscous  perturbation  can  reach 
the  same  order  of  magnitude  as  the  potential  flow  velocity;  i.e., 
0(U  *)~0(Ua  ) .  Hence,  this  method  can  treat  medium  separation. 

“  O  “0 

The  solution  thus  obtained,  including  the  circulation  zone,  is 
rational.  In  the  two-dimensional  case,  the  results  agree  with 
other  computational  results  and  experimental  data. 

As  compared  to  conventional  methods  to  calculate 


viscous/inviscid  interaction,  the  time  and  capacity  required  by 
this  method  are  significantly  less  because  there  is  no  need  to 
repeat  the  inviscid  portion  in  the  iterations.  The  use  of 
integral  boundary  layer  equation  also  helps  to  reduce  the  demand 
on  the  computer.  We  believe  that  it  can  be  applied  to  a  variety 
of  complex  flow  situation  due  to  its  simplicity  and  adaptability. 

Appendix:  Based  on  reference  [10],  the  two-dimensional 
incompressible  interacting  boundary  layer  equations  are: 


(Ai) 


dX 

dx 


dH 

dX 


dd*  f.d6  ,fidH  dX 
~d^  =  Hd7+e-dTd7 
dU,_  p  .  dd* 

~d T-R^+a"~dT 


where 


^-=4t41  -  c"  *  ,  C„  - - l- - - 

2  U  ,0  *■ 


R.  =  pR,*  ,  -  + — 4* - 


2  2 


A.  =  OJ,j 5*);,  (U,t6*)s  , 


r  -~U-^-+dR'  dd<'-  y  a  ^ 

*'-^r+-dT  £  ,:d.x 


xir 


h  is  the  step.  (  )q  and'  (  )N  represent  the  starting  point 

and  the  ending  point,  respectively. 
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Numerical  Computation  of  Extended  Kalman  Filter  and  Its 


Application  to  Aerodynamic  Parameter  Identification  of  a  Reentry 


Satellite 


Chen  Qiongkang  and  Jiang  Quanwei 


(China  Aerodynamic  Research  and  Development  Center) 


Abstract 


This  paper  describes  several  numerical  computation  problems 


to  identify  aerodynamic  parameters  by  using  an  extended  Kalman 


filter.  The  ultimate  goal  is  to  save  computer  time  and  memory  to 


the  extent  possible  and  to  improve  the  speed  of  computation.  It 


is  discussed  in  the  following  parts:  optimization  of  numerical 


algorithm  to  solve  the  covariance  matrix  differential  equations. 


computation  of  the  Kalman  gain  and  its  correction,  an  algorithm 


to  calculate  the  process  noise  in  order  to  save  computer  time, 


and  an  independent  verification  method  using  numerical 


derivatives.  These  methods  have  been  applied  to  the 


identification  of  aerodynamic  parameters  based  on  the  reentry 


flight  data  of  a  certain  satellite  and  the  results  obtained  are 


satisfactory. 


1.  Introduction 


To  use  an  extended  Kalman  filter  to  identify  aerodynamic 


parameters  by  numerical  computation,  there  are  generally  three 


requirements: 


a# 


■VVVV4".  • 


(1)  Fast  Computing  Speed.  Because  a  great  deal  of 
computation  is  required  for  model  selection,  performance  study 
and  error  analysis,  increasing  computing  speed  has  significance 
in  both  computing  efficiency  and  economic  benefit.  If  we 
investigate  real  time  estimation  of  aerodynamic  properties  by  an 
on-board  computer,  then  computing  speed  and  memory  storage  become 
especially  important.  In  this  paper,  an  optimized  numerical 
algorithm  to  solve  the  covariance  differential  matrix  equation 
and  to  calculate  the  Kalman  gain  and  its  correction  is  introduced 
in  Sections  2  and  3.  As  compared  to  unoptimized  algorithms,  the 
speed  can  be  improved  by  five  fold  and  the  storage  can  be  reduced 
by  three  fourths. 

(2)  Good  Filter  Characteristics.  Due  to  various  errors 
such  as  mathematical  model,  measurement  and  numerical 
computation,  the  aerodynamic  parameter  identified  by  an  extended 
Kalman  filter  is  divergent  in  some  cases.  In  addition,  if  there 
is  no  process  noise,  the  filter  is  more  sensitive  to 
computational  error.  In  the  fourth  section  of  this  paper  a 
method  to  add  process  noise  is  introduced.  The  focus  is  placed 
on  a  numerical  method  to  calculate  this  process  noise  to  save 
computing  time.  It  is  a  variable  separation  method. 

(3)  Accuracy  of  Filter  Computation.  In  the  extended  Kalman 
filter  computation  for  a  large  non-linear  system  with  several 
variables,  it  is  necessary  to  calculate  a  large  number  of  partial 
derivatives.  Their  derivations  are  more  complicated  and  mistakes 
are  easily  made.  Next,  the  quality  of  the  filter,  to  a  large 
extent,  depends  on  the  variations  of  these  partial  derivatives  in 


their  variable  range  (dynamic  variables  and  parameters) .  In  the 
fifth  section  of  this  paper,  a  numerical  method  is  presented  to 
calculate  the  partial  derivatives  which  can  independently  verify 
the  accuracy  of  an  analytical  partial  derivative.  It  also  can 
reflect  the  linearity  of  the  dynamic  system. 

In  this  work,  we  combined  our  effort  reported  in  reference 
[1]  and  referred  to  recent  foreign  publications  [2],  [3],  [4], 

[5]  and  [6]  to  study  the  problems  mentioned  above.  Our  optimized 
numerical  method  is  explained.  A  FORTRAN  program  was  compiled. 

Many  results  were  obtained  on  a  Model  320  computer.  Moreover, 

Manuscript  received  on  February  25,  1984,  revision  received  on 
September  16. 

it  has  already  been  applied  to  aerodynamic  parameter  / 9 7 

identification  of  satellite  reentry  with  satisfactory  results. 

2.  Optimized  Numerical  Method  for  Covariance  Matrix 
Differential  Equations 

The  non-linear  random  differential  equation  to  determine  the 
state  of  the  system  is  given  as  follows; 

X(t)  *  f IX ( t) ,  t]  +  W ( t)  (1) 

The  measurement  equation  is  a  non-linear  function  of  the  state 
vector  X(t)  in  a  discrete  form; 

-K  =  — K ^ ^  +  -K 

*10 


(2) 


where  W(t)  is  the  Gauss  white  noise  whose  mean  value  is  zero. 


the  frequency  spectrum  density  is  Q(t).  VR  is  the  independent 


Gauss  random  sequence  whose  mean  value  is  zero.  The  covariance 


matrix  is  RK  and  K  is  the  sampling  sequence.  The  state  vector 


X ( t )  comprised  of  the  dynamic  quantity  Y(t)  and  the  extended 


coefficient  parameter  C(t)  ,  i.e. 


Therefore,  equation  (1)  can  also  be  written  as 


*(')==  Q  (0 


£<0|  f£CK(O.C,0|  f  Wit) 


where  g[Y(t),  C,  t]  is  a  non-linear  vector  function.  The 


superscript  represents  a  derivative  with  respect  to  time  t. 


The  subscript  represents  a  vector. 


In  extended  Kalman  filter  computation,  it  is  necessary  to 


solve  a  non-linear  ordinary  differential  equation  and  an  error 


covariance  matrix  differential  equation: 


X(t)  =  f  [X ( t)  ,  t] 


P(t)  =  F  ( t)  P  ( t)  +  P(t)FT(t)  +  Q  ( t) 


where 


r-  bf- 

F(/)  —  IV  A 

OA  ,!-*«> 


-  I-  ' 


which  is  called  the  state  matrix.  P(t)  is  a  non-negative  fixed 


matrix.  The  superscript  "T"  represents  a  matrix  transformation 


X(t)  is  the  estimated  value  of  the  state  vector. 


In  order  to  save  computing  time  and  storage,  two  important 


mmm 


properties  of  equation  (6),  i.e.,  the  sparseness  of  the  state 
matrix.  F(t)  and  the  symmetry  of  P(t)  and  P(t),  are  used  to 
optimize  the  numerical  method  to  solve  equation  (6) .  It  is 
illustrated  as  follows: 

The  state  matrix  F(t)  is  a  sparse  matrix.  Its  non-zero 
elements  are  distributed  in  1-S  rows;  i.e.: 

F(i)  =  \Fss  F*‘l  (8) 

1  0  0  J 

Here,  Fsg(S  x  an<^  FS  x  a)  are  sParse  matrices  where  S  + 
a=N.  S  and  a  are  the  dimensionals  of  Y  and  C,  respectively.  N 
is  the  order  of  the  state  matrix  F(t).  First,  we  employed  the 
method  reported  in  reference  [3]  to  solve  equation  (6)  by  blocks. 
It  takes  less  time  than  other  conventional  methods.  However,  the 
objective  of  maximum  savings  in  computing  time  and  storage  is  not 
met.  To  this  end,  we  used  the  characteristics  of  the  sparse 
matrix.  The  algorithm  is  to  compress  M  non-zero  elements  into  a 
series  of  one-dimensional  numbers  AN.  Furthermore,  two  series  IA 
and  JA  are  used  to  describe  the  data  structure  of  the  matrix 
F ( t) .  AN  and  JA  have  M  components  each.  I A  has  N  +  1 
components.  Hence,  there  are  2M  +  N  +  1  storage  elements  in  this 
algorithm^7! .  Based  on  the  relation  between  these  components  in 
matrix  operation,  it  is  possible  to  find  F(t)*P(t).  By  treating 
the  non-zero  elements  of  the  matrix  F(t)  as  vectors,  the  number 
of  operations  is  the  least. 

Based  on  the  calculated  F(t)P(t)  and  the  symmetry  of  P(t) , 
we  get  P(t)FT(t)  =  [F(t) P(t) ]T,  i.e.,  P(t)FT(t)  can  be  determined 

Aik. 


by  transformation.  With  the  process  noise  Q(t) ,  and  numerical 
integration,  we  can  find  the  solution  P(t)  of  the  covariance 
matrix  equation  (6) .  Due  to  the  symmetry  of  equation  (6) ,  all 
calculations  were  made  in  the  upper  triangle  of  the  matrix. 

In  the  meantime,  equation  (5)  is  also  numerically 
integrated.  The  integrated  between  two  sampling  points  must  be 
done  in  several  steps.  The  solutions  of  equation  (5)  and  (6)  at 
the  Kfch  sampling  point  are  denoted  as  XK(-)  and  PK(-), 

respectively  (K=l,2, . ),  which  are  used  as  the  estimated 

values. 

3.  Kalman  Gain  and  Its  Correction 

The  estimated  state  quantity  and  covariance  matrix  XR(-)  and 

A 

PK(~)  must  be  corrected  to  obtain  the  corrections  XK(+)  and  PR(+) 
at  the  K*-*1  sampling  point.  The  equations  are: 

kk  =  pk(-)hJ[hkpk(-)hI+rk1“1  O) 

XK(  +  )  =  X^  ( -  )  +Kj(  [  Zpr-hj^ ( X^  ( -  )  )  ]  (10) 

PK(+)  =  (I-KkHk)Pk(-)  (11) 

where  KR  is  the  Kalman  gain  matrix.  I  is  an  order  unity 
matrix. 


dht(X) 
Hl  =  ~~dX~ 


I  * 

jr-?*'-’ 


(12) 


It  is  a  measurement  matrix,  as  well  as  a  sparse  matrix. 


This  section  is  concentrated  with  matrix  operations.  The 
same  method  as  used  in  the  previous  section  is  used  to  store  and 


The 


compute  the  non-zero  elements  of  the  sparse  matrix  HK. 
matrix  is  obtained  and  stored.  Because  equation  (9)  is 

used  in  two  places,  its  transpose  will  be  used  in  equation  (11) . 

Due  to  the  symmetry  of  PR(-)  and  RK,  computing  [HKPK(-) 
hJ+Rk]  is  carried  out  on  one  side  of  the  main  diagonal  line.  Its 
inverse  is  obtained  by  Cholesky  decomposition  in  this  paper.  For 
a  positive  definite  matrix  A,  we  have  A=LLT.  Here,  L  is  a  lower 
triangular  matrix.  Based  on  A“^= (LLT) “*= (L-1) TL_1 ,  after 
decomposing  the  matrix  A,  it  is  only  required  to  find  the  inverse 
of  the  lower  triangular  matrix.  The  amount  of  computation  is 
relatively  small. 

In  order  to  simplify  the  program  and  to  reduce  the  number  of 
cycles,  after  the  inverse  matrix  [HRPK (-) h^+Rk] -1  is  obtained,  we 
can  sequentially  calculate  equations  (9) ,  (10)  and  (11)  by  using 
the  results  of  each  column  vector  of  KR.  A  FORTRAN  program  was 
written  which  consists  of  one  external  loop  with  four  small 
loops. 


4.  A  Time-saving  Process  Noise  Computation  Method 

This  work  distributes  six  process  noises  into  three  force 
equations  and  three  moment  equations.  The  error  of  this 
mathematical  model  will  be  numerically  integrated  with  other 
equations.  Together  with  measurement  noise,  the  simulated 
observation  data  for  filter  identification  are  obtained.  Here, 
the  process  noise  is  considered  as  a  stable  Gauss  process  which 
is  simulated  by  the  following  method^  : 
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*  • 

irr(o=  2  ^“(©sJ+^.k) 

*• » 


(13) 


i  =  I? 2f . .  6 

W^(t)  represents  the  process  noise,  a'^K~N(0 ,crfj  )  r  v ir  is  a  random 
number  evenly  distributed  over  0~2tp,  ©e=(K-1/2)  x  ©r/Nq,  ©£#  N0, 
are  known. 

It  is  obvious  that  equation  (13)  is  the  sum  of  NQ  terms  of 
the  sine  series  of  the  integral  variable  t.  If  we  computer  W^(t) 
directly,  it  takes  a  great  deal  of  work  in  the  forward  problem. 

In  order  to  save  computer  time,  this  article  introduces  a 
variable  separation  method  to  calculate  the  process  noise  (13) . 
Specifically,  equation  (13)  is  expanded  to  separate  the 
trigonometric  function  terms  containing  the  time  variable  t  from 
those  that  do  not  contain  t;  i.e. 

A’  « 

W ,  (/)  =  £)(a, c  cos  q),t  sin  ©xf  +  o,*  sin  ip, t  cos  ©*<> 

i- 1  •  ( 14) 

i  -  1,  2  . ,6 

After  the  expansion,  C^Kcos<^K,  <a^Ksin«P£K  are  independent  of 
time,  and  can  be  calculated  and  stored.  When  numerically 
integrating  equation  (1)  ,  sinw^t,  cosco^t  are  two  trigonometric 
functions  related  to  time  t.  But,  they  are  independent  of  the 
function  number  i.  The  number  of  times  to  call  for  sine  and 
cosine  functions  is  one  third  of  that  when  directly  calculating 
equation  (13).  In  other  words,  this  variable  separation  method 
can  save  two  thirds  of  the  computer.  . 


5.  Verification  by  Numerical  Derivatives 


As  described  above,  two  partial  derivative  matrices  (7)  and 
(12)  must  be  calculated  in  filter  equations.  The  partial 
derivatives  are  calculated  using  an  analytic  method.  Its 
derivation  is  more  complicated  and  it  is  more  mistake  prone. 

In  order  to  verify  the  accuracy  of  the  analytical  partial 
derivative,  and  to  demonstrate  the  linearity  between  the  non¬ 
linear  functions  f  and  h,  numerical  partial  derivative  is  sought 
in  this  work.  For  example,  with  respect  to  the  matrix  F,  the 
computation  formula  is: 

„  /.(Jr  +  AA’,>-/.(jf-AA’()  ,,r, 

F  =  “  iKjrr  (15) 

When  AXj  is  chosen  to  be  sufficiently  small,  equation  (15)  can 
accurately  approximate  the  analytical  partial  derivative.  When 
AXj  is  chosen  to  be  the  initial  standard  deviation  of  the  state 
quantity,  it  is  possible  to  show  the  linearity  of  the  function  f 
in  the  neighborhood  of  the  initial  standard  deviation. 

Therefore,  we  choose  AXj  =  a  \/  (Pj j)  o  where  a  is  a  proportionality 
factor  which  is  between  0 . 01~1 .0,  \/(Pj j ) q  is  the  initial  standard 
deviation. 

The  computation  shows  that  when  a  =  0.01~0.1  the  analytical 
partial  derivative  agrees  with  the  numerical  partial  derivative. 
When  a  =  1.0,  the  results  are  very  close.  This  shows  that  the 
function  f  has  a  good  linearity  in  the  neighborhood  of  the  chosen 
standard  deviation.  Similar  computation  can  be  made  for  Hi^. 


6.  Conclusions 


This  numerical  method  has  been  used  in  the  identification  of 
aerodynamic  parameters  based  on  the  reentry  flight  data  of  a 
certain  satellite.  Major  aerodynamic  properties  of  the  stars 
were  identified  based  on  high  and  low  altitude  data, 
respectively,  including  axial  force,  static  stability  and  non- 
symmetric  moment  coefficient.  The  results  are  satisfactory. 

In  this  work,  we  referenced  the  technical  report  by  Comrade 
Xu  Jinzhi  on  conventional  methods.  We  wish  to  express  our 
gratitude. 
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Relaxation  Phenomenon  Behind  A  Bow  Shock  Wave  In  Dusty  Gas  / 1 0 1 
Zhao  Guoying  (Institute  of  Mechanics,  Chinese  Academy  of  Sciences) 
Zhing  Xichang  (Center  of  Computation,  Chinese  Academy  of  Sciences) 

The  supersonic  flow  behind  the  explosion  wave  often  contains 
dust.  When  a  Pitot  pressure  gage  is  used  to  measure  its  total 
pressure,  the  effect  of  the  dusty  motion  should  be  considered. 

When  an  object  is  traveling  at  high  speed  in  an  atmosphere 
containing  raindrops,  ice  particles  and  fog,  the  effect  of  solid 
and  liquid  micro  particles  should  also  be  taken  into  account. 

This  type  of  problem  is  called  the  supersonic  dusty  flow.  Its 
flow  pattern  is  as  follows:  A  bow  shock  wave  is  formed  in  front 
of  the  dull  body.  The  flow  velocity  drops  to  zero  when  reaching 
the  surface.  When  a  micro  particle  passes  through  the  shock  wave 
plane,  its  parameters  remain  unchanged.  Afterward,  due  to  its 
retardation  by  the  gas  and  the  heat  exchange  between  them,  the 
particle  velocity  decreases  and  its  temperature  rises.  When  it 
reaches  the  surface,  the  dust  particle  still  has  certain  normal 
velocity.  They  will  bounce  off  the  surface  and  be  smashed.  They 
can  even  erode  the  surface. 

When  analyzing  this  type  of  flow,  most  authors  neglected  the 
effect  of  the  particles  on  the  gas  dynamics  and  used  the 
trajectory  theory  to  perform  numerical  analysis.  Only  Changt^ 
considered  the  effect  of  the  particles  on  the  flow.  However,  the 
direct  method  and  series  expansion  technique  used  by  him  is  only 
applicable  to  high  Mach  number  flows.  In  addition,  because  the 
drag  coefficient  equation  is  wrong,  his  results  are  only 


qualitative.  Besides,  his  temperature  relaxation  process  is  not 
even  qualitatively  accurate. 

In  this  work,  a  two-way  flow  model  and  the  straight  line 
method  is  used  to  study  the  relaxation  phenomenon  numerically. 

The  relaxation  process  of  a  dusty  gas  behind  a  blow  shock  wave  is 
discussed  in  detail. 


1.  Basic  Equations  and  Boundary  Conditions 


Let  us  assume  that  the  gas  is  an  ideal  gas  whose  isobaric 
heat  capacity  Cp  and  specific  heat  Y  are  constants.  The  dust 
particles  are  identical  spheres,  evenly  distributed  in  the  flow. 

In  addition,  the  inside  of  the  sphere  is  at  a  uniform  temperature. 
The  size  of  the  dust  particle  is  far  larger  than  the  mean  free 
path  of  the  gas  molecules.  Therefore,  compared  to  the  momentum 
transfer  between  the  dust  particles  and  the  gas,  the  shock  wave 
thickness  of  the  gas  can  be  neglected.  The  particles  occupy  a 
very  small  percentage  of  the  volume.  Collision  between  particles 
can  be  neglected.  The  amount  of  particles  is  sufficiently  high  in 
the  volume  of  interest.  Consequently,  it  can  be  treated  as  a 
continuous  medium.  There  are  only  momentum  and  energy  exchanges 
between  the  dust  and  the  gas.  There  is  no  mass  transfer. 

Under  the  above  assumptions,  the  axisymmetric  flow  equation 
of  a  dusty  gas  without  inviscid  stress  effect  can  be  written  as 
follows  by  using  the  boundary  layer  coordinate  (See  Figure  1) : 


(A) 
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where  o  is  the  gas  concentration.  hx=l+y/Rj3  where  Rb  is  the 
radius  of  curvature  of  the  surface  and  R  is  the  gas  constant. 


the  dust 


( otu,r )  -f  —^(.o,vprh.)  *=0 


-F„ 

-T„ 


For 


Here,  ap  is  the  dust  particle  concentration.  The  mutual 
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where  T^d^  Ps/3MCDRe,  Re=aV^(up-u)  2+  ( vp-v)  2  dp/p. 

where  dp  is  the  particle  diameter  and  ps  is  its  density.  As  for 

the  drag  coefficient  CD,  the  spherical  drag  coefficient  standard 


curve  is  used. 


-Ir-  (1  +  0.0625/?e  — 3.508  X  10" '/?e:-r  8.914  x  10":i?e'  ),  Re<  180 
Re 


-|*-(3.1  +  0.0167/?e>  , 
Re 


Re>  180 


The  heat  exchange  between  the  two  planes  is 


Qt=otC,(T,-T)lTT 


where  Cs  is  the  specific  heat  of  the  particle  material. 
TT  -  PsCs/6kNu 

where  k  is  the  thermal  conductivity  of  the  gas. 

Nu  =  2.0  +  0.459Re°*55Pr0,33 


The  work  done  by  the  dust  particles  on  the  gas  is 
wk  55  °p[(Up-U)2  +  (Vp-v)2]/^ 


k  -  wpL  V  Up-u^  ■  w  m 

The  boundary  conditions  for  the  above  equations  are:  the 


parameters  behind  the  bow  shock  wave  are  given  by  Hugoniot 


correlations. 


h  +  Vll2-hm  +  Vl,s/2  , 


p+pV*  =  pm  +  cmVi,N 


“"--('-•ir)’’-''-"  • 


i  mm 


where  the  subscript  N  is  the  normal  direction  of  the  shock  wave 
plane.  nx  and  ny  are  the  components  of  the  unit  normal  vector  of 
the  wave  front  in  x  and  y  directions,  respectively.  The  flow  / 1 03 

condition  for  the  solid  particle  behind  the  bow  shock  wave  is: 

Up=U®r  Vp=Vm 

On  the  surface,  the  normal  velocity  of  the  gas  is  zero?  i.e., 
vJy_Q=o.  As  for  the  solid  particles,  let  us  assume  that  they  are 
absorbed  by  the  surface  after  they  collide  with  it.  We  do  not 
consider  their  burning,  crushing  and  eroding  processes. 

2.  Solving  the  Equations,  Results  and  Discussion 

Let  us  assume  that  the  object  is  a  sphere  and  the  gas  is  air. 
Let  us  use  the  radius  of  the  sphere  Rb,  incoming  flow  velocity  .Vcor 
and  \^/R  as  the  characteristic  scales  to  make  length, 
velocity,  concentration  (density) ,  pressure  and  temperature 
dimensionless,  respectively*.  Let  us  assume  that  the  shock  wave 
surface  equation  is  y=f^(x)  and  the  surface  equation  is  y=f2(x)=0. 
Let  us  introduce  the  following  transformations  £=y/fi(x)  and  *)=x. 

We  solved  the  above  boundary  value  problem  by  a  linear  method. 

The  major  concept,  details  and  problems  of  the  linear  method  are 
discussed  in  reference  [2] . 


Figure  1.  Shock  Wave  Shapes  and  Sonic  Speed  Streamlines  at 
Various  Particle  Concentrations 

1.  shock  wave  plane 

2.  sonic  speed  streamline 

3.  surface 


Figure  2.  Pressure  Distribution  on  Surface  vs.  Particle 
Concentration 


In  the  stationary  point  zone,  although  a  source  term  is  added 
to  the  right  side  of  the  equation  of  motion  due  to  the  dust 
effect,  the  equation  is  still  hybrid  in  nature.  Therefore,  there 
is  no  difficulty  to  use  the  straight  line  method  to  convert  it 
into  a  boundary  value  problem  of  a  regular  differential  equation. 
The  dust  particle  equation  is  parabolic.  Thus,  as  long  as  we  do 
not  consider  particles  bouncing  off  the  surface,  there  is  no 
problem  to  integrate  it  by  the  straight  line  method. 

In  this  work,  the  above  equations  were  solved  by  using  the 
straight  line  method  over  a  wide  range  of  Mach  number,  McO(2~10^)  , 
and  Opco(0~1.0).  When  there  is  order  of  magnitude  difference 
between  Bj (=3CD/4dpps)  and  B2  (“SNu/Re^PrCp) ,  due  to  the  large 
difference  between  the  relaxation  distance  of  temperature  and  that 
of  velocity,  a  stiff  problem  occurs  in  these  equations.  In  our 
program,  we  employed  the  Treanor  method to  eliminate  this 
problem.  We  choose  the  situation  B^=4. 5,  B2=1.0. 

(corresponding  to  Re00=6932,  NUoc=54.6),  Op^* 0—1.0,  Cp/Cg=1.0  and  Pr 
=0.69,  the  numerical  results  are  shown  in  Figures  1-7  as  well  as 
in  the  table.  Figure  1  shows  the  shock  wave  shapes  and  sonic 
streamlines  at  various  CTp^.  The  table  lists  the  variation  of 
separation  distance  A  with  As  increases,  the  shock  wave 

continues  to  get  close  to  the  surface.  The  sonic  streamline  moves 
upstream  near  the  surface.  When 

♦In  the  text  below,  with  the  expectation  of  quantities  with  *, 
others  are  dimensionless. 


dp®  is  small,  the  effect  of  particle  motion  on  the  gas  flow  field  / 1 04- 
is  not  significant.  Figure  2  shows  the  variation  of  stationary 
point  and  surface  pressure  at  different  o_  .  when  increases, 

k*CO  pcD 

the  surface  pressure  also  increases.  When  ^  =i.o,  the  surface 

pressure  increases  one  fold.  Figure  3  shows  the  normal  velocity 
of  the  gas  and  that  of  the  particle.  Although  the  relaxation 
process  between  them  is  similar  to  that  in  the  normal  shock  wave 
case^4! ,  and  the  only  difference  is  located  far  away  from  the 
stationary  streamline,  despite  both  the  gas  and  the  particle  having 
transversal  velocity,  yet  the  gas  flow  field  has  a  transversal 
pressure  gradient  while  it  is  not  true  for  the  particle  flow 
field.  Therefore,  the  relaxation  process  is  much  more  complicated 
as  compared  to  the  normal  shock  wave  case.  Figure  4  shows  the 
transversal  velocity  relaxation  process  using  ®pffl“1.0  as  an 
example.  Far  away  from  the  stationary  streamline,  the  transversal 
velocity  relaxation  distance  is  shorter  than  that  near  the 
stationary  point.  This  is  due  to  the  effect  of  the  gas  pressure 
gradient.  Near  the  surface,  the  transversal  velocity  of  the  gas 
is  higher  than  that  of  the  particle.  The  particle  is  carried  by 
the  gas. 


Table  Shock  Wave  Separation  Distance  vs.  cr*  in  Dusty  Air 
When  M  3-  =  3  p®  . 


Or. 
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Figure  3.  Variation  of  Normal  Velocities  of  Gas  and  Particle 


1.  normal  velocity 

2.  gas  velocity 

3.  particle  velocity 


Figure  4.  Variation  of  Transversal  Velocities  of  Gas  and 
Particle 

1.  transversal  velocity 


Figure  5.  Gas  and  Particle  Temperature  Variation  vs.  o 

1.  dimensionless  temperature 

2.  gas  temperature  T 

3.  particle  temperature  T 


Figure  5  shows  the  relaxation  process  between  the  temperature 
fields  of  these  two  phases.  Due  to  the  fact  that  this  is  a  slow 
process  and  the  heat  transferred  from  the  gas  to  the  particle 
behind  the  shock  wave  front  is  less  than  the  frictional  work  done 
by  the  particle  on  the  gas,  the  temperature  rises  when  crpm  is  not 
equal  to  zero.  The  particle  temperature,  however,  does  not  vary 
significantly.  (In  the  figure,  the  particle  temperature  is  low 
when  cTp^l.O,  is  attributed  to  the  thinning  of  the  shock  wave 
layer.)  This  is  different  from  Chang's  results.  Figure  6  shows  / 105 
the  particle  concentration  variation  along  the  stationary 
streamline.  It  is  apparent  that  its  value  is  increased  five  fold. 
Away  from  the  stationary  streamline,  the  concentration  increase  is 
also  more  or  less  similar.  Figure  7  shows  the  effect  of  the  Mach 
number  Mm  and  on  the  separation  distance  under  the  condition 
of  B2  and  B2  as  described  before. 


Figure  6.  Variation  of  Particle  Concentration  Along  Stationary 
Streamline  at  Various  o_, 

p  CO 

1.  particle  concentration 


Figure  7.  Variation  of  Separation  Distance  A  With  Different  M( 


a.*f» 


Our  computation  shows  that  the  relaxation  process  between 
these  two  phases  is  related  to  and  B2.  B^  controls  the  dynami 
relaxation  process  and  B2  affects  temperature  relaxation.  In  the 
example  given  in  this  work,  if  B^  and  B2  are  increased,  the 
relaxation  process  is  gradually  shifted  near  the  shock  wave.  If 
B^  and  B2  are  decreased,  the  coupling  of  these  two  phases  is 
reduced.  As  far  as  the  effect  of  the  particle  on  the  gas  flow 
field  is  concerned,  we  can  use  <7_  Bi  and  cr  as  two  combined 

parameters.  Only  when  they  are  very  small,  we  can  solve  the  flow 
field  in  the  stationary  zone  by  single  coupling.  The  entire 
computation  shows  that  as  long  as  B^  and  B2  are  equal,  and 
remains  the  same,  the  variation  of  the  relaxation  processes 
plotted  in  Figures  2~6  versus  Mm  is  very  small. 

The  method  introduced  here  can  be  expanded  to  non-spher ical 
head  objects  and  non-axisymmetric  flow.  Non-equilibrium  effects 
such  as  chemical  reactions  in  the  gas  phase  and  phase 
transformation  of  the  particles  can  also  be  concluded. 

The  authors  wish  to  express  their  gratitude  to  comrade  Yu 
Hongru  for  his  guidance. 
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A  Simple  Formula  of  Drag  Coefficient  of  Sphere  in  Rarefied  / 1 0 7 

Hypersonic  Flow 
Ma  Jiahuan 

(Institute  of  Mechanics,  Academia  Sinica) 

Abstract 

This  paper  presents  a  simple  relationship  of  the  spherical 
drag  coefficient  using  the  viscous  parameter  K=lk/V  ReDj  a 
variable  which  is  suited  for  rarefied  hypersonic  flow.  In  the  ReD 
=  5  x  102~105  and  =  8~25  range,  it  agrees  with  the  experimental 
data.  In  addition,  the  dynamic  pressures  of  calibrating  spheres 
desired  from  the  spherical  drag  coefficient  based  on  this  formula 
agree  with  those  obtained  by  Pitot  pressure  measurements. 

1.  Introduction 

The  drag  coefficient  of  a  sphere  in  a  rarefied  hypersonic 
flow  is  a  problem  of  great  concern  because  it  has  many  practical 
applications  such  as  in  estimating  the  drag  and  life  of  satellites 
and  in  high  level  atmospheric  density  measurements.  We  are 
concerned  about  this  issue  because  the  dynamic  pressure  is  usually 
obtained  by  using  a  collaborating  sphere  when  a  model  is  tested  in 
a  hypersonic  pulse  wind  tunnel.  Thus,  the  drag  coefficient  of  the 
sphere  under  the  experimental  conditions  must  be  known  ahead  of 
time.  Hence,  it  is  necessary  to  provide  a  simple  and  accurate 
spherical  drag  coefficient  formula. 

2R 


As  we  all  know,  in  a  continuous  hypersonic  flow,  the 
spherical  drag  coefficient  is  almost  a  constant.  However,  the 
flow  region  is  not  the  concern.  Once  the  flow  enters  the  so- 
called  skidding  and  transition  regions,  the  spherical  drag 
coefficient  not  only  varies  with  the  Reynolds  number  but  also 
relates  to  the  Mach  number.  In  these  flow  regions,  many 
experiments  have  been  conducted  abroadl®-^.  Nevertheless,  at 
high  Mach  numbers,  these  formulas  cannot  totally  reflect  the 
variation  of  the  drag  coefficient.  To  this  end,  this  paper 
presents  a  simple  formula  using  the  viscosity  parameter  K  =  M^ReD 
as  a  variable. 

2.  Presentation  of  the  Formula 

* 

The  flow  regions  are  divided  based  on  the  viewpoint  expressed 
in  reference  [12] .  The  mean  free  molecular  path  to  boundary  layer 
thickness  ratio  is  used  to  express  the  magnitude  of  viscosity. 
Regions  with  different  flow  characteristics  are  thus  defined  as 
follows.  When  the  mean  free  molecular  path  to  boundary  layer 
thickness  ratio  is  x/a<0.01,  it  is  a  continuous  flow.  When  0.01< 
x/gd,  it  is  in  a  skidding  region.  When  X/s=-0,  it  is  a  free 
molecular  flow.  Between  the  skidding  zone  and  free  molecular 
flow,  it  is  a  transition  zone.  These  regions  are  divided  as  shown 
in  Figurel.  In  addition,  we  also  included  the  various  model 
experimental  states  tested  in  a  coordinate  system  comprised  of  the 
Reynolds  number  and  Mach  number  of  a  free  flow  using  the 
calibrating  sphere  diameter  as  a  characteristic  length.  In  most 


cases,  as  -we  can  see,  the  free  flying  calibrating  sphere  does  not 
belong  to  the  continuous  flow  category.  Of  course,  it  is  not 
proper  if  the  same  spherical  drag  coefficient  is  selected  as  in 
the  continuous  flow  case. 

Manuscript  received  on  March  12,  1984.  Revision  received  on 
September  29. 

Table  1.  Comparison  of  Dynamic  Pressures  Obtained  by  Free  Flying  / 1 08 
Calibrating  Sphere  to  Pitot  Measurements 


1.  equipment 

2.  qs(kg/cm“)  obtained  from  calibrating  sphere 
3*  q»(kg/cm1 2)  obtained  from  Pitot  pressure 


Figure  1.  Division  of  Flow  Regions 

1.  free  molecular  flow 

2.  transition  region 

3.  skidding  region 

4.  continuous  flow 

How  should  we  choose  the  spherical  drag  coefficient  value  in 
this  region?  Based  on  the  recommendation  of  reference  [10] ,  we 
made  corrections  by  an  increment  ACD«1/ yReD  in  the  skidding  zone 
(ReD=102~104) .  A  great  deal  of  experimental  data  was  summarized 
in  that  paper.  However,  most  experimental  results  were  obtained 
in  supersonic  range.  The  formula  of  the  drag  coefficient  of  a 
sphere  is  given  as  follows: 

CD  =  0.95  +  5/-/Re^  (1) 

where  ReD  is  the  Reynolds  number  of  the  flow  based  on  a 
characteristic  length  equal  to  the  diameter  of  the  sphere  and  CD 
is  the  total  drag  coefficient  of  the  sphere.  The  linking  of  the 
experimental  value  of  the  drag  coefficient  in  hypersonic 
continuous  flow  is  considered  primarily  based  on  the  result 
reported  in  reference  [13]  =  CD=0.915.  Hence,  we  recommended  the 
following  formula  in  reference  [14], 

CD  -  0.899  +  5/  y^ReD 


(2) 


It  is  obvious  that  this  equation  cannot  reflect  the  effect  of  Mach 
number.  When  the  Reynolds  number  is  still  relatively  large,  such 
as  ReDe>10*,  the  effect  of  viscosity  is  a  small  correction  term 
because  inertia  is  still  the  primary  force;  the  formula  can  still 
be  used  without  causing  major  deviations.  However,  in  a 
hypersonic  flow,  when  the  Reynolds  number  is  even  smaller,  the 
limitation  of  this  equation,  same  as  that  of  equation  (1),  is  even 
more  obvious. 

Let  us  introduce  the  spherical  drag  relationship  presented  in 
reference  [11]  which  is  also  a  result  of  an  analysis  of  a  great 
deal  of  experimental  data,  particularly  the  target  data  measured 
by  A.B.  Bailey  and  J.  Hiatt it  is  a  new  relationship  derived 
based  on  a  theoretical  analysis  to  calculate  the  drag  coefficient 
of  a  sphere.  It  incorporates  the  effect  of  wall  temperature  and 
considers  the  applicability  of  the  entire  flow  region.  In  the 
hypersonic  region,  it  can  be  expressed  as: 
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Where  S  is  the  molecular  speed  ratio;  =  M^,.  Vy  /2;  V  is  the 
specific  heat  of  the  gas;  Tw  is  the  wall  temperature  of  the  sphere 
and  Tro  is  the  temperature  of  the  incident  flow. 

This  formula  takes  various  factors  influencing  the  drag 


coefficient  of  the  sphere  into  account.  Nevertheless,  its 
applicable  range  is  limited  to  the  supersonic  region  where  M®<6. 
In  the  hypersonic  case,  let  us  first  examine  the  variation  of  the 


spherical  drag  coefficient  in  the  region  (See  Figure  2) .  This 
figure  comes  from  reference  [9]  which  uses  ReD  as  a  variable 
parameter.  It  shows  the  pattern  how  drag  coefficient  varies  Mach 
number.  The  solid  lines  are  the  distribution  patterns  obtained  by 
a  large  number  of  data  points. 


Figure  2.  Hypersonic  Spherical  Drag  Coefficient  Over  a  Large 
Range 

1.  free  molecular  flow  limit 

As  we  can  see,  when  Mc0<  6,  the  drag  coefficient  of  the  sphere 
decreases  as  increases.  This  is  consistent  with  the  trend 
given  by  equation  (3).  After  M^IO,  however,  the  drag  coefficient 
increases  as  increases.  This  trend  is  especially  obvious  when 
ReD  is  very  small.  In  this  case,  equation  (3)  cannot  reflect  this 
trend.  We  believe  that  the  reason  why  the  drag  coefficient  of  a 
sphere  in  a  rarefied  hypersonic  flow  increases  rapidly  with 


increasing  Mach  number  is  because  the  shock  wave  is  adhered  closer 
to  the  surface  at  hypersonic  speeds.  In  addition,  because  the  gas 
stream  is  rarefied,  the  boundary  layer  thickness  is  larger. 

Therefore,  a  viscous  interference  between  the  shock  wave  and  the 
boundary  layer,  similar  to  that  frequently  found  in  the  hypersonic 
flow  around  a  dull  slender  body,  will  take  place.  Consequently, 
the  drag  coefficient  is  significantly  increased.  For  this  reason, 
we  suggest  to  use  the  viscous  parameter  K=Maj/-y/ ReD  to  express  the 
drag  coefficient,  similar  to  the  division  of  flow  characteristic 
region.  Through  an  order  of  magnitude  analysis,  in  the  range  of 
concern  ReD>>l,  we  know  that  the  viscous  parameter  K  is 
proportional  to  X/g.  We  use  K=M*,/ \ZRej)  to  express  the  variation  / 110 
of  drag  coefficient  of  a  sphere  in  a  rarefied  hypersonic  flow. 
Furthermore,  we  consider  the  fact  that  there  should  be  no 
discontinuities  in  this  formula  in  the  case  of  continuous  flow 
where  the  viscosity  effect  is  minimal.  Thus,  the  following 
equation  is  given: 


■0.9 


.  1 


M. 


2  N/tf 


(4) 


eD 


The  applicable  range  of  this  formula  is 

M.  =  8~25 ,  ReD  =  5  x  lO^-lO5  (5) 

It  corresponds  to  K=0.025~1.10.  The  formula  is  simple  and  neat. 

It  can  reflect  the  variation  of  drag  coefficient  in  hypersonic 
rarefied  flow  over  a  large  range.  The  trend  of  increase  in  the  M. 
drag  coefficient  calculated  based  on  this  formula  agrees  with  that 
reflected  by  a  great  deal  of  experiments  at  high  Mach  numbers. 


Only  when  ReD  =  1  x  103,  its  deviation  becomes  significant.  The 
maximum  deviation  can  reach  14%.  In  other  situations,  the 
deviation  is  much  smaller.  Especially  when  ReD>104,  the 
calculated  results  agree  well  with  the  experimental  points. 

Now,  let  us  compare  the  spherical  drag  coefficients  in 
hypersonic  rarefied  flows  measured  in  shock  tube  wind  tunnels  in 
references  [7]  and  [15]  to  the  values  calculated  based  on  this 
equation.  They  are  shown  in  Figure  3  using  the  viscous  parameter 
K  =  M*,/  VReD  as  the  coordinate.  The  straight  line  in  the  figure  .3 
is  derived  from  equation  (4)  which  is  very  close  to  the 
experimental  values  obtained  at  four  Mach  numbers. 


Figure  3.  Comparison  of  Spherical  Drag  Coefficient  to 
Experimental  Data 

1.  equation  (4)  in  this  paper 


In  addition,  the  drag  coefficients  derived  from  equation  (4) 
also  agree  with  the  dynamic  pressures  obtained  with  free  flying 
calibrating  spheres  and  Pitot  force  measurements  in  shock  tubes 
and  cannon  wind  tunnels  (See  Table  1) . 

Finally,  it  should  be  pointed  out  that  the  ratio  of  the  wall 
temperature  of  the  sphere  to  the  flow  temperature  has  an  effect  on 
its  drag  coefficient.  Reference  [16]  was  dedicated  to  studying 
this  problem.  Its  result  shows  that  the  effect  of  wall 
temperature  is  relatively  small  when  the  Reynolds  number  is  high. 
In  the  applicable  range  of  the  equation,  as  long  as  the  wall 
temperature  is  not  much  higher  than  the  flow  temperature,  this 
effect  can  be  neglected.  As  for  other  factors,  such  as  the  effect 
of  sphere  material  and  surface  condition  investigated  in  reference 
[7],  no  influence  was  found  within  the  accuracy  of  measurement.  . 
As  a  simple  formula,  the  effect  of  secondary  factors  is  not 
reflected  in  this  work. 

3.  Conclusions 

In  this  work  a  simple  formula  to  calculate  the  drag 
coefficient  of  a  sphere  in  hypersonic  rarefied  flow  is  introduced. 
A  viscous  parameter,  k  =  M^/  ^ ReD,  is  used  as  a  convection  factor 
to  reflect  the  changes  of  the  drag  coefficient  with  Re  and  Hw.  In 
the  applicable  range  of  the  formula:  Ks  =  8~25;  ReD  =  5  x  10^~ 
10^,  the  calculated  results  are  in  agreement  with  a  vast  number  of 
experimental  data. 

Furthermore,  the  formula  also  agrees  with  the  dynamic 
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pressures  measured  by  calibrating  spheres  and  with  Pitot 
pressures.  We  believe  that  this  formula  offers  a  convenient  way 
to  calculate  the  drag  coefficient  with  good  accuracy  when  we  make 
measurement  of  free  flights  in  rarefied  hypersonic  flows. 
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Alleviation  and  Control  of  Asymmetric  Load  at  High  Angle  of  Attack 

Yang  Yongnian 

(Northwestern  Polytechnical  University) 

Abstract 

This  paper  introduces  several  methods  to  alleviate  and 
control  asymmetric  load  at  high  angles  cf  attack,  as  well  as  their 
research  status  and  typical  results.  Comments  on  active  control 
of  asymmetric  load  in  a  wind  tunnel  are  presented. 

Symbols 

‘a  angle  of  attack  Cy  lateral  force  coefficient 

aonset  onset  angle  Cn  yawing  moment  coefficient 

Ci  rolling  moment  coefficient  C„  blowing  coefficient 

1.  Introduction 

In  recent  years  a  great  deal  of  work  on  the  aerodynamic 
characteristics  at  high  attack  angles  has  been  done  by  scholars  in 
the  field.  In  particular,  a  large  number  of  reports  have  been 
published  to  describe  the  occurrence  and  development  of 
-„.r^otr^c  ioa<3  high  angles  of  attack,  and  the  visualization  and 
measurement  of  flow  patterns  in  leeward  areas,  through 
experimental  research.  The  results  of  these  studies  show  that 
when  an  aircraft  is  flying  symmetrically  at  a  high  angle  of 
attack,  asymmetric  loads  which  vary  according  to  a  complicated 
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pattern  will  occur  which  is  unfavorable.  In  addition,  after 
stalling,  the  rudder  efficiency  drops,  the  maneuverability  and 
operability  of  the  aircraft  will  deteriorate.  In  this  case,  the 
asymmetric  load  at  high  attack  angles  may  be  used  to  improve  the 
lateral  maneuverability  and  controllability  of  the  aircraft. 

Thus,  the  subject  of  alleviation  and  control  of  asymmetric  load 
was  brought  up.  It  is  necesary  to  find  an  appropriate  way  not 
only  to  alleviate  the  deleterious  effect  of  asymmetric  load  at 
high  angles  of  attack  but  also  to  improve  the  performance  of  the 
aircraft.  Various  methods  of  alleviation  and  control  were 
presented.  The  concept  of  active  control  was  introduced.  Due  to 
the  complexity  of  this  problem,  it  is  still  in  the  feasibility 
stage.  Nevertheless,  it  is  a  topic  of  importance  for  further 
research. 

This  paper  reviews  several  methods  presented  in  recent  years 
to  alleviate  and  control  asymmetric  load  at  high  angle  of  attack 
and  presents  their  typical  results.  These  methods  are: 
installing  transition  strips  at  the  nose,  installing  rings  or 
fixed  strakes  at  the  nose,  changing  nose  shape,  rotating  nose, 
controlled  strakes  and  blowing  vortex  control.  Finally,  some 
ideas  on  active  control  in  wind  tunnels  are  also  presented. 

The  methods  mentioned  above  and  the  typical  results  are 
briefly  described  as  follows: 

2.  Transition  Strip 

The  most  commonly  used  method  to  alleviate  the  asymmetric 


load  at  high  angles  of  attack  is  to  install  a  transition  strip  at 
a  suitable  place  on  the  nose  to  change  the  attached  layer  to  a 
turbulent  attached  layer  to  reduce  the  separation  zone.  Thus,  the 
onset  of  the  asymmetric  vortex  in  the  leeward  region  is  delayed 
and  its  asymmetry  is  weakened .  Typical  results  are  shown  in 
Figure  la. 

Manuscript  received  on  June  14,  1984,  revision  received  on 
September  18. 

3.  Ring  on  the  Nose  or  Fixed  Strake  / 113 

Reference  [5]  reported  the  wind  tunnel  experimental  results 
of  the  nose  of  a  rotating  body  installed  with  a  ring  or  a  fixed 
strake.  The  results  showed  that  asymmetric  load  could  be  reduced 
by  adding  a  ring  or  a  fixed  strake  at  the  nose  and  the  onset  angle 
was  increased.  The  experimental  results  are  shown  in  Figures  lb 
and  1c,  respectively. 
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Figure  1.  Typical  Experimental  Results  of  Transition  Strip,  Ring 
and  Fixed  Strake 


1.  without  transition  strip 

2.  with  transition  strip 

3.  without  ring 

4.  with  ring 

5.  without  strake 

6.  with  strake 


4.  Changing  Nose  Shape 


References  [6]  and  [7]  pointed  out  that  the  geometric  shape 


of  the  nose  has  an  important  effect  on  the  asymmetric  load.  Major 


parameters  affecting  the  asymmetric  load  include:  aspect  ratio  of 


the  nose,  baseline  shape,  radius  of  curvature  at  the  tip  of  the 


nose,  and  the  cross-sectional  shape  of  the  nose.  Experimental 


results  show  that  reducing  the  aspect  ratio  of  the  nose, 


air 


increasing  the  radius  of  curvature  at  the  tip,  increasing  the 


concave  curvature  of  the  baseline  and  increasing  the  horizontal 


width  of  the  nose  (i.e.  locating  the  major  axis  of  the  ellipse 


horizontally)  can  reduce  the  asymmetric  load  and  increase  the 


onset  angle.  Typical  results  are  shown  in  Figure  2. 
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Figure  2.  Effect  of  Nose  Shape  on  Asymmetric  Load 


1.  radius  of  nose  tip 

2.  base  radius 

3.  aspect  ratio 

4.  cone 

5.  normal  arch 

6.  aircraft 


5.  Rotating  Nose  Device 


Reference  [3]  introduced  the  concept  of  rotating  nose.  A 
portion  of  the  nose  is  rotating  at  a  constant  speed  in  a  fixed 
direction  around  its  axis  to  reduce  or  even  eliminate  the 
asymmetry  in  leeward  flow.  Experimental  results  pointed  out  that 
the  asymmetric  lateral  force  could  be  significantly  reduced  when 
the  rotating  speed  is  higher  than  8  rev/sec.  A  schematic  diagram 
of  this  rotating  device  and  its  typical  results  are  shown  in 
Figure  3.  Its  effect  on  the  pitch  moment  is  not  very  high. 

6.  Controlled  Strake 

Reference  [8]  introduced  the  concept  of  hinged  strake.  The 
strake  is  connected  to  the  body  by  hinges,  with  varying  angle  of 
attack,  the  strake  angle  can  vary  symmetrically  or  asymmetrically 
according  to  a  specific  pattern  to  alleviate  the  deleterious 
effect  caused  by  the  breakdown  of  strake  vortex  to  improve  the 
aerodynamic  properties  of  the  strake.  The  results  showed  that  it 
is  possible  to  find  a  better  reverse  angle  on  the  strake  (-30°  on 
both  right  and  left  strakes)  so  that  the  lateral  force  and  yawing 
moment  are  essentially  zero  in  the  attack  angle  range  of  0^-40° 

(See  Figure  4a) .  In  addition,  the  longitudinal  characteristics 
are  also  improved.  The  experimental  results  show  that  asymmetric  /1 14 
yawing  of  the  strakes  can  provide  higher  rolling  moments  (See 
Figure  4b) .  Therefore,  using  .asymmetric  strakes  is  an  effective 
device  for  lateral  force  control  of  aircraft  at  high  angle  of 


attack 


Figure  3.  Schematic  Diagram  of  Rotating  Nose  and  Its  Typical 
Results 

1.  rotating  part 
2..  motor 

3.  forward 

4.  reverse 

5.  forward 

6.  reverse 
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Figure  4.  Typical  Results  of  Asymmetric  Strakes 


1.  without  strake 

2.  without  strake 


7.  Active  Blowing  Device 


References  [4],  [9],  [10]  and  [11]  presented  various  schemes 


of  blowing  in  various  directions  (along  the  surface  tangentially 


downstream  or  upstream,  blowing  in  normal  direction)  to  destroy 


the  asymmetric  flow  in  the  leeward  region  to  make  it  symmetric. 


Thus,  the  asymmetric  load  is  alleviated  or  eliminated.  Or,  its 


direction  even  changed.  Figures  5a  and  b  show  the  typical  results 


r 
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of  normal  blowing^1 2 3 4 5 6 7!  and  tangential  blowing^10!,  respectively. 

The  results  show  that  blowing  is  very  effective  in  alleviating 
asymmetric  load.  Reference  [11]  pointed  out  that  the  axial 
position  of  the  nozzle  should  be  placed  as  close  to  the  nose  as 
possible.  The  circumferential  position  of  the  nozzle  should  be  at 
the  onset  of  the  separation  vortex.  Thus,  the  blowing  efficiency 
can  be  improved. 


Figure  5.  Typical  Blowing  Results 

1.  without  blowing 

2.  blowing  upstream 

3.  blowing  downstream 

4.  normal  direction 

5.  tangential  direction 

6.  angle  of  attack 

7.  semi-conical  angle 


8.  Comments 

In  the  methods  discussed  above ,  the  alleviation  of  asymmetric 

load  by  the  nose  transition  strip,  ring  and  fixed  strake,  and 

changing  geometric  shape  of  the  nose  cannot  be  automatically 

adjusted  based  on  the  flight  condition.  Therefore,  they  belong  to 

the  passive  control  category.  These  methods  are  relatively  simple 

and  effective  in  specific  ranges.  More  studies  are  currently 

being  conducted  on  these  methods  as  well.  Rotating  nose, 

controlled  strake  and  active  blowing  device  can  automatically  be  / 115 

changed  based  on  different  flight  conditions.  They  belong  to  the 

active  control  category.  These  methods  are  more  complex  and  are 

more  difficult  to  realize.  They  are  still  in  a  research  stage.  A 

* 

great  deal  of  theoretical  analyses  and  experimental  work  will  be 
needed.  The  following  specific  problems  must  be  addressed  in 
performing  experiments  on  active  control  systems  at  high  attack 
angles  in  wind  tunnels: 

1.  We  must  understand  the  correlation  between  the  active 
control  parameters  and  the  asymmetric  load  in  order  to  optimize 
the  parameters  of  the  active  control  device. 

2.  In  the  design  and  adjustment  of  the  control  system,  we 
must  resolve  problems  such  as  the  selection  of  sensors  and 
execution  mechanism,  determining  the  control  law,  and  the  design, 
fabrication  and  adjustment  of  the  control  circuit. 

3.  In  terms  of  model  design,  because  sensors  and  active 
control  device  must  be  installed  on  the  model,  model  design  is 
very  difficult. 


4.  In-  formulating  the  experimental  protocol,  we  should  first 


consider  the  control  with  a  single  factor,  test  its  feasibility 
and  resolve  all  the  technical  problems.  Then,  we  will  consider 
the  overall  effect  on  the  performance  of  the  aircraft.  In 
addition,  after  an  active  control  device  is  added,  the  movement  of 
the  load  and  the  control  device  becomes  a  dynamic  problem.  It  is 
required  to  consider  the  dynamic  effect  of  the  model.  In  this 
case,  rigid  degree  of  freedom  of  the  aircraft  must  be  simulated. 

In  conclusion,  active  control  of  the  asymmetric  load  at  high 
anqles  of  attack  is  a  technically  complicated  but  rewarding 
research  subject.  It  is  still  in  its  infancy.  It  should  attract 
people's  attention  for  further  exploration. 
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